Nek5000
SEM for Incompressible NS
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
fast3d.F90
Go to the documentation of this file.
1 !-----------------------------------------------------------------------
3 subroutine gen_fast_spacing()
4  use input, only : param
5  use speclib, only : zwgll, zwgl
6  implicit none
7 
8 #if 0
9  use kinds, only : dp
10  real(DP) :: x(nx1,ny1,nz1,nelv)
11  real(DP) :: y(nx1,ny1,nz1,nelv)
12  real(DP) :: z(nx1,ny1,nz1,nelv)
13 
14  use size_m
15  use parallel
16  use soln
17  use wz_m
18 
19  integer, parameter(lxx=lx1*lx1)
20 
21  common /ctmpf/ lr(2*lx1+4),ls(2*lx1+4),lt(2*lx1+4) &
22  , llr(lelt),lls(lelt),llt(lelt) &
23  , lmr(lelt),lms(lelt),lmt(lelt) &
24  , lrr(lelt),lrs(lelt),lrt(lelt)
25  real(DP) :: lr ,ls ,lt
26  real(DP) :: llr,lls,llt
27  real(DP) :: lmr,lms,lmt
28  real(DP) :: lrr,lrs,lrt
29 
30  integer :: lbr,rbr,lbs,rbs,lbt,rbt,e
31 
32  real(DP) :: axwt(lx2)
33 
34  ierr = 0
35 #endif
36  if (param(44) == 1) then
37  write(*,*) "Oops: param(44)"
38 #if 0
39  ! __ __ __
40  ! Now, for each element, compute lr,ls,lt between specified planes
41 
42  n1 = nx2
43  n2 = nx2+1
44  nz0 = 1
45  nzn = 1
46  if (if3d) then
47  nz0= 0
48  nzn=n2
49  endif
50  eps = 1.e-7
51  if (wdsize == 8) eps = 1.e-14
52 
53  ! Find mean spacing between "left-most" planes
54  call plane_space2(llr,lls,llt, 0,wxm2,x,y,z,n1,n2,nz0,nzn)
55 
56  ! Find mean spacing between "middle" planes
57  call plane_space(lmr,lms,lmt, 1,n1,wxm2,x,y,z,n1,n2,nz0,nzn)
58 
59  ! Find mean spacing between "right-most" planes
60  call plane_space2(lrr,lrs,lrt,n2,wxm2,x,y,z,n1,n2,nz0,nzn)
61 #endif
62  else
63  call load_semhat_weighted ! Fills the SEMHAT arrays
64  endif
65 
66  return
67 end subroutine gen_fast_spacing
68 
69 !-----------------------------------------------------------------------
71 subroutine plane_space(lr,ls,lt,i1,i2,w,x,y,z,nx,nxn,nz0,nzn)
72  use kinds, only : dp
73  use size_m, only : nelv
74  use input, only : if3d
75  implicit none
76 
77  integer :: nx, nz0, nxn, nzn, i1, i2
78  real(DP) :: w(*),lr(*),ls(*),lt(*)
79  real(DP) :: x(0:nxn,0:nxn,nz0:nzn,*)
80  real(DP) :: y(0:nxn,0:nxn,nz0:nzn,*)
81  real(DP) :: z(0:nxn,0:nxn,nz0:nzn,*)
82  real(DP) :: lr2,ls2,lt2
83 
84  integer :: ny, nz, j1, k1, j2, k2, ie, k, j, i
85  real(DP) :: wsum, weight
86 
87 ! Now, for each element, compute lr,ls,lt between specified planes
88  ny = nx
89  nz = nx
90  j1 = i1
91  k1 = i1
92  j2 = i2
93  k2 = i2
94 
95  do ie=1,nelv
96 
97  if (if3d) then
98  lr2 = 0.
99  wsum = 0.
100  do k=1,nz
101  do j=1,ny
102  weight = w(j)*w(k)
103  ! lr2 = lr2 + ( (x(i2,j,k,ie)-x(i1,j,k,ie))**2
104  ! $ + (y(i2,j,k,ie)-y(i1,j,k,ie))**2
105  ! $ + (z(i2,j,k,ie)-z(i1,j,k,ie))**2 )
106  ! $ * weight
107  lr2 = lr2 + weight / &
108  ( (x(i2,j,k,ie)-x(i1,j,k,ie))**2 &
109  + (y(i2,j,k,ie)-y(i1,j,k,ie))**2 &
110  + (z(i2,j,k,ie)-z(i1,j,k,ie))**2 )
111  wsum = wsum + weight
112  enddo
113  enddo
114  lr2 = lr2/wsum
115  lr(ie) = 1./sqrt(lr2)
116 
117  ls2 = 0.
118  wsum = 0.
119  do k=1,nz
120  do i=1,nx
121  weight = w(i)*w(k)
122  ! ls2 = ls2 + ( (x(i,j2,k,ie)-x(i,j1,k,ie))**2
123  ! $ + (y(i,j2,k,ie)-y(i,j1,k,ie))**2
124  ! $ + (z(i,j2,k,ie)-z(i,j1,k,ie))**2 )
125  ! $ * weight
126  ls2 = ls2 + weight / &
127  ( (x(i,j2,k,ie)-x(i,j1,k,ie))**2 &
128  + (y(i,j2,k,ie)-y(i,j1,k,ie))**2 &
129  + (z(i,j2,k,ie)-z(i,j1,k,ie))**2 )
130  wsum = wsum + weight
131  enddo
132  enddo
133  ls2 = ls2/wsum
134  ls(ie) = 1./sqrt(ls2)
135 
136  lt2 = 0.
137  wsum = 0.
138  do j=1,ny
139  do i=1,nx
140  weight = w(i)*w(j)
141  ! lt2 = lt2 + ( (x(i,j,k2,ie)-x(i,j,k1,ie))**2
142  ! $ + (y(i,j,k2,ie)-y(i,j,k1,ie))**2
143  ! $ + (z(i,j,k2,ie)-z(i,j,k1,ie))**2 )
144  ! $ * weight
145  lt2 = lt2 + weight / &
146  ( (x(i,j,k2,ie)-x(i,j,k1,ie))**2 &
147  + (y(i,j,k2,ie)-y(i,j,k1,ie))**2 &
148  + (z(i,j,k2,ie)-z(i,j,k1,ie))**2 )
149  wsum = wsum + weight
150  enddo
151  enddo
152  lt2 = lt2/wsum
153  lt(ie) = 1./sqrt(lt2)
154 
155  else ! 2D
156  lr2 = 0.
157  wsum = 0.
158  do j=1,ny
159  weight = w(j)
160  ! lr2 = lr2 + ( (x(i2,j,1,ie)-x(i1,j,1,ie))**2
161  ! $ + (y(i2,j,1,ie)-y(i1,j,1,ie))**2 )
162  ! $ * weight
163  lr2 = lr2 + weight / &
164  ( (x(i2,j,1,ie)-x(i1,j,1,ie))**2 &
165  + (y(i2,j,1,ie)-y(i1,j,1,ie))**2 )
166  wsum = wsum + weight
167  enddo
168  lr2 = lr2/wsum
169  lr(ie) = 1./sqrt(lr2)
170 
171  ls2 = 0.
172  wsum = 0.
173  do i=1,nx
174  weight = w(i)
175  ! ls2 = ls2 + ( (x(i,j2,1,ie)-x(i,j1,1,ie))**2
176  ! $ + (y(i,j2,1,ie)-y(i,j1,1,ie))**2 )
177  ! $ * weight
178  ls2 = ls2 + weight / &
179  ( (x(i,j2,1,ie)-x(i,j1,1,ie))**2 &
180  + (y(i,j2,1,ie)-y(i,j1,1,ie))**2 )
181  wsum = wsum + weight
182  enddo
183  ls2 = ls2/wsum
184  ls(ie) = 1./sqrt(ls2)
185  ! write(6,*) 'lrls',ie,lr(ie),ls(ie)
186  endif
187  enddo
188  return
189 end subroutine plane_space
190 
191 !-----------------------------------------------------------------------
192 subroutine get_fast_bc(lbr,rbr,lbs,rbs,lbt,rbt,e,bsym,ierr)
193  use size_m, only : ndim
194  use input, only : cbc
195  use parallel, only : lglel
196  use topol, only : eface
197  use tstep, only : ifield
198  implicit none
199 
200  integer :: lbr,rbr,lbs,rbs,lbt,rbt,e,bsym
201  integer :: fbc(6)
202  integer :: iface, ied, ibc, ierr
203 
204 ! ibc = 0 <==> Dirichlet
205 ! ibc = 1 <==> Dirichlet, outflow (no extension)
206 ! ibc = 2 <==> Neumann,
207 
208  do iface=1,2*ndim
209  ied = eface(iface)
210  ibc = -1
211 
212 !max if (ifmhd) call mhd_bc_dn(ibc,iface,e) ! can be overwritten by 'mvn'
213 
214  if (cbc(ied,e,ifield) == ' ') ibc = 0
215  if (cbc(ied,e,ifield) == 'E ') ibc = 0
216  if (cbc(ied,e,ifield) == 'msi') ibc = 0
217  if (cbc(ied,e,ifield) == 'MSI') ibc = 0
218  if (cbc(ied,e,ifield) == 'P ') ibc = 0
219  if (cbc(ied,e,ifield) == 'p ') ibc = 0
220  if (cbc(ied,e,ifield) == 'O ') ibc = 1
221  if (cbc(ied,e,ifield) == 'ON ') ibc = 1
222  if (cbc(ied,e,ifield) == 'o ') ibc = 1
223  if (cbc(ied,e,ifield) == 'on ') ibc = 1
224  if (cbc(ied,e,ifield) == 'MS ') ibc = 1
225  if (cbc(ied,e,ifield) == 'ms ') ibc = 1
226  if (cbc(ied,e,ifield) == 'MM ') ibc = 1
227  if (cbc(ied,e,ifield) == 'mm ') ibc = 1
228  if (cbc(ied,e,ifield) == 'mv ') ibc = 2
229  if (cbc(ied,e,ifield) == 'mvn') ibc = 2
230  if (cbc(ied,e,ifield) == 'v ') ibc = 2
231  if (cbc(ied,e,ifield) == 'V ') ibc = 2
232  if (cbc(ied,e,ifield) == 'W ') ibc = 2
233  if (cbc(ied,e,ifield) == 'SYM') ibc = bsym
234  if (cbc(ied,e,ifield) == 'SL ') ibc = 2
235  if (cbc(ied,e,ifield) == 'sl ') ibc = 2
236  if (cbc(ied,e,ifield) == 'SHL') ibc = 2
237  if (cbc(ied,e,ifield) == 'shl') ibc = 2
238  if (cbc(ied,e,ifield) == 'A ') ibc = 2
239  if (cbc(ied,e,ifield) == 'S ') ibc = 2
240  if (cbc(ied,e,ifield) == 's ') ibc = 2
241  if (cbc(ied,e,ifield) == 'J ') ibc = 0
242  if (cbc(ied,e,ifield) == 'SP ') ibc = 0
243 
244  fbc(iface) = ibc
245 
246  if (ierr == -1) write(6,1) ibc,ied,e,ifield,cbc(ied,e,ifield)
247  1 format(2i3,i8,i3,2x,a3,' get_fast_bc_error')
248 
249  enddo
250 
251  if (ierr == -1) call exitti('Error A get_fast_bc$',e)
252 
253  lbr = fbc(1)
254  rbr = fbc(2)
255  lbs = fbc(3)
256  rbs = fbc(4)
257  lbt = fbc(5)
258  rbt = fbc(6)
259 
260  ierr = 0
261  if (ibc < 0) ierr = lglel(e)
262 
263 ! write(6,6) e,lbr,rbr,lbs,rbs,(cbc(k,e,ifield),k=1,4)
264 ! 6 format(i5,2x,4i3,3x,4(1x,a3),' get_fast_bc')
265 
266  return
267  end subroutine get_fast_bc
268 
269 !-----------------------------------------------------------------------
274 subroutine load_semhat_weighted() ! Fills the SEMHAT arrays
275  use size_m, only : nx1
276  use semhat, only : ah, bh, ch, dh, zh, dph, jph, bgl, zgl, dgl, jgl, wh
277  implicit none
278  integer :: nr
279 
280  nr = nx1-1
281  call generate_semhat(ah,bh,ch,dh,zh,dph,jph,bgl,zgl,dgl,jgl,nr,wh)
282  call do_semhat_weight(jgl,dgl,bgl,nr)
283 
284  return
285 end subroutine load_semhat_weighted
286 
287 !----------------------------------------------------------------------
289 subroutine do_semhat_weight(jgl,dgl,bgl,n)
290  use kinds, only : dp
291  implicit none
292  integer :: n
293  real(DP) :: bgl(1:n-1),jgl(1:n-1,0:n),dgl(1:n-1,0:n)
294  integer :: i, j
295  do j=0,n
296  do i=1,n-1
297  jgl(i,j)=bgl(i)*jgl(i,j)
298  enddo
299  enddo
300  do j=0,n
301  do i=1,n-1
302  dgl(i,j)=bgl(i)*dgl(i,j)
303  enddo
304  enddo
305  return
306 end subroutine do_semhat_weight
307 
308 !-----------------------------------------------------------------------
328 subroutine generate_semhat(a,b,c,d,z,dgll,jgll,bgl,zgl,dgl,jgl,n,w)
329  use kinds, only : dp
330  use speclib, only : zwgll, zwgl
331  implicit none
332 
333  integer :: n
334  real(DP) :: a(0:n,0:n),b(0:n),c(0:n,0:n),d(0:n,0:n),z(0:n)
335  real(DP) :: dgll(0:n,1:n-1),jgll(0:n,1:n-1)
336 
337  real(DP) :: bgl(1:n-1),zgl(1:n-1)
338  real(DP) :: dgl(1:n-1,0:n),jgl(1:n-1,0:n)
339 
340  real(DP) :: w(0:(n+1)*2)
341 
342  integer :: np, nm, n2
343  integer :: i, j, k
344 
345  np = n+1
346  nm = n-1
347  n2 = n-2
348 
349  call zwgll(z,b,np)
350 
351  do i=0,n
352  call fd_weights_full(z(i),z,n,1,w)
353  do j=0,n
354  d(i,j) = w(j+np) ! Derivative matrix
355  enddo
356  enddo
357 
358  if (n == 1) return ! No interpolation for n=1
359 
360  do i=0,n
361  call fd_weights_full(z(i),z(1),n2,1,w(1))
362  do j=1,nm
363  jgll(i,j) = w(j ) ! Interpolation matrix
364  dgll(i,j) = w(j+nm) ! Derivative matrix
365  enddo
366  enddo
367 
368  a = 0._dp
369  do j=0,n
370  do i=0,n
371  do k=0,n
372  a(i,j) = a(i,j) + d(k,i)*b(k)*d(k,j)
373  enddo
374  c(i,j) = b(i)*d(i,j)
375  enddo
376  enddo
377 
378  call zwgl(zgl,bgl,nm)
379 
380  do i=1,n-1
381  call fd_weights_full(zgl(i),z,n,1,w)
382  do j=0,n
383  jgl(i,j) = w(j ) ! Interpolation matrix
384  dgl(i,j) = w(j+np) ! Derivative matrix
385  enddo
386  enddo
387 
388  return
389 end subroutine generate_semhat
390 
391 !-----------------------------------------------------------------------
409 subroutine fd_weights_full(xx,x,n,m,c)
410  use kinds, only : dp
411  implicit none
412 
413  integer :: n, m
414  real(DP) :: xx, x(0:n),c(0:n,0:m)
415 
416  integer :: i, j, k, mn
417  real(DP) :: c1, c2, c3, c4, c5
418 
419  c1 = 1.
420  c4 = x(0) - xx
421 
422  do k=0,m
423  do j=0,n
424  c(j,k) = 0.
425  enddo
426  enddo
427  c(0,0) = 1.
428 
429  do i=1,n
430  mn = min(i,m)
431  c2 = 1.
432  c5 = c4
433  c4 = x(i)-xx
434  do j=0,i-1
435  c3 = x(i)-x(j)
436  c2 = c2*c3
437  do k=mn,1,-1
438  c(i,k) = c1*(k*c(i-1,k-1)-c5*c(i-1,k))/c2
439  enddo
440  c(i,0) = -c1*c5*c(i-1,0)/c2
441  do k=mn,1,-1
442  c(j,k) = (c4*c(j,k)-k*c(j,k-1))/c3
443  enddo
444  c(j,0) = c4*c(j,0)/c3
445  enddo
446  c1 = c2
447  enddo
448 ! call outmat(c,n+1,m+1,'fdw',n)
449  return
450 end subroutine fd_weights_full
451 
452 !-----------------------------------------------------------------------
453 subroutine swap_lengths()
454  use kinds, only : dp
455  use size_m, only : lx1, ly1, lz1, lelv
456  use size_m, only : nx1, nelv
457  use geom, only : xm1, ym1, zm1
458  use hsmg, only : llr, lrr, lmr, lls, lms, lrs, llt, lmt, lrt
459  use input, only : if3d
460  use wz_m, only : wxm1
461  implicit none
462 
463  real(DP), allocatable :: l(:,:,:,:)
464  integer :: e, n2, nz0, nzn, nx, n, j, k
465 
466  allocate(l(lx1, ly1, lz1, lelv))
467 
468  n2 = nx1-1
469  nz0 = 1
470  nzn = 1
471  nx = nx1-2
472  if (if3d) then
473  nz0 = 0
474  nzn = n2
475  endif
476  call plane_space(lmr,lms,lmt,0,n2,wxm1,xm1,ym1,zm1,nx,n2,nz0,nzn)
477 
478  n=n2+1
479  if (if3d) then
480  do e=1,nelv
481  do j=2,n2
482  do k=2,n2
483  l(1,k,j,e) = lmr(e)
484  l(n,k,j,e) = lmr(e)
485  l(k,1,j,e) = lms(e)
486  l(k,n,j,e) = lms(e)
487  l(k,j,1,e) = lmt(e)
488  l(k,j,n,e) = lmt(e)
489  enddo
490  enddo
491  enddo
492  call dssum(l)
493  do e=1,nelv
494  llr(e) = l(1,2,2,e)-lmr(e)
495  lrr(e) = l(n,2,2,e)-lmr(e)
496  lls(e) = l(2,1,2,e)-lms(e)
497  lrs(e) = l(2,n,2,e)-lms(e)
498  llt(e) = l(2,2,1,e)-lmt(e)
499  lrt(e) = l(2,2,n,e)-lmt(e)
500  enddo
501  else
502  do e=1,nelv
503  do j=2,n2
504  l(1,j,1,e) = lmr(e)
505  l(n,j,1,e) = lmr(e)
506  l(j,1,1,e) = lms(e)
507  l(j,n,1,e) = lms(e)
508  ! call outmat(l(1,1,1,e),n,n,' L ',e)
509  enddo
510  enddo
511  ! call outmat(l(1,1,1,25),n,n,' L ',25)
512  call dssum(l)
513  ! call outmat(l(1,1,1,25),n,n,' L ',25)
514  do e=1,nelv
515  ! call outmat(l(1,1,1,e),n,n,' L ',e)
516  llr(e) = l(1,2,1,e)-lmr(e)
517  lrr(e) = l(n,2,1,e)-lmr(e)
518  lls(e) = l(2,1,1,e)-lms(e)
519  lrs(e) = l(2,n,1,e)-lms(e)
520  enddo
521  endif
522  return
523 end subroutine swap_lengths
524 
525 !----------------------------------------------------------------------
527 subroutine row_zero(a,m,n,e)
528  use kinds, only : dp
529  implicit none
530 
531  integer :: m,n,e
532  real(DP) :: a(m,n)
533  integer :: j
534  do j=1,n
535  a(e,j)=0.
536  enddo
537  return
538 end subroutine row_zero
539 !-----------------------------------------------------------------------
541 subroutine swap(b,ind,n,temp)
542  use kinds, only : dp
543  implicit none
544  real(DP) :: B(*),TEMP(*)
545  integer :: n, IND(*)
546  integer :: i, jj
547 
548 !***
549 !*** SORT ASSOCIATED ELEMENTS BY PUTTING ITEM(JJ)
550 !*** INTO ITEM(I), WHERE JJ=IND(I).
551 !***
552  DO i=1,n
553  jj=ind(i)
554  temp(i)=b(jj)
555  END DO
556  DO i=1,n
557  b(i)=temp(i)
558  END DO
559  RETURN
560 end subroutine swap
561 
row_zero
subroutine row_zero(a, m, n, e)
zero the eth row of a
Definition: fast3d.F90:527
dssum
subroutine dssum(u)
Direct stiffness sum.
Definition: dssum.F90:54
get_fast_bc
subroutine get_fast_bc(lbr, rbr, lbs, rbs, lbt, rbt, e, bsym, ierr)
Definition: fast3d.F90:192
tstep
cleaned
Definition: tstep_mod.F90:2
input
Input parameters from preprocessors.
Definition: input_mod.F90:11
soln
Definition: soln_mod.F90:1
wz_m
cleaned
Definition: wz_mod.F90:2
gen_fast_spacing
subroutine gen_fast_spacing()
Generate fast diagonalization matrices for each element.
Definition: fast3d.F90:3
hsmg
Module containing data for HSMG.
Definition: hsmg_mod.F90:3
do_semhat_weight
subroutine do_semhat_weight(jgl, dgl, bgl, n)
re-weight jgl and dgl with bgl
Definition: fast3d.F90:289
topol
Arrays for direct stiffness summation. cleaned.
Definition: topol_mod.F90:3
semhat
cleaned
Definition: semhat_mod.F90:2
fd_weights_full
subroutine fd_weights_full(xx, x, n, m, c)
This routine evaluates the derivative based on all points in the stencils. It is more memory efficien...
Definition: fast3d.F90:409
parallel::lglel
integer function lglel(iel)
Definition: parallel_mod.F90:176
swap_lengths
subroutine swap_lengths()
Definition: fast3d.F90:453
load_semhat_weighted
subroutine load_semhat_weighted()
Note that this routine performs the following matrix multiplies after getting the matrices back from ...
Definition: fast3d.F90:274
generate_semhat
subroutine generate_semhat(a, b, c, d, z, dgll, jgll, bgl, zgl, dgl, jgl, n, w)
Generate matrices for single element, 1D operators:
Definition: fast3d.F90:328
parallel
cleaned
Definition: parallel_mod.F90:2
geom
Geometry arrays.
Definition: geom_mod.F90:2
speclib::zwgll
subroutine, public zwgll(Z, W, NP)
Definition: speclib.F90:95
swap
subroutine swap(b, ind, n, temp)
Reorder vector using temporary buffer.
Definition: fast3d.F90:541
plane_space
subroutine plane_space(lr, ls, lt, i1, i2, w, x, y, z, nx, nxn, nz0, nzn)
Here, spacing is based on harmonic mean. pff 2/10/07.
Definition: fast3d.F90:71
exitti
subroutine exitti(stringi, idata)
Definition: comm_mpi.F90:328
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
Generated on Tue Jul 14 2015 10:52:58 for Nek5000 by   doxygen 1.8.8