Nek5000
SEM for Incompressible NS
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
bdry.F90
Go to the documentation of this file.
1 !-----------------------------------------------------------------------
3 SUBROUTINE setlog()
4  use size_m, only : ndim, nelv, nelt, nfield, nid, lxd
5  use ctimer, only : ifsync
6  use geom, only : ifvcor, ifgeom, ifsurt, ifwcno, ifeppm, ifqinp, ifmelt
7  use input, only : param, ifadvc, iftmsh, ifneknek, ifmoab, ifkeps, ifmodel
8  use input, only : ifnonl, cbc, ifheat, ifmvbd, ifflow, ifnav, ifstrs, ifprint
9  use input, only : ifaxis, ifcyclic, ifchar, ifusermv, ifuservp, iflomach
10  use input, only : ifsplit, iftran, nobj, lochis, ifintq, hcode, nhis
11  use tstep, only : ifield, nelfld
12  use turbo, only : ifswall, ifcwuz
13  implicit none
14 
15  CHARACTER(3) CB
16  LOGICAL :: IFALGN,IFNORX,IFNORY,IFNORZ
17 
18  integer :: nface, nmxv, nmxt, iel, ifc, iq, ih, iobj
19 
20  nface = 2*ndim
21  nmxv = nface*nelv
22  nmxt = nface*nelt
23 
24  ifprint = .true.
25  ifvcor = .true.
26  ifgeom = .false.
27  ifintq = .false.
28  ifsurt = .false.
29  ifwcno = .false.
30  ifswall = .false.
31  DO ifield=1,nfield
32  ifnonl(ifield) = .false.
33  END DO
34 
35  CALL lfalse(ifeppm,nmxv)
36  CALL lfalse(ifqinp,nmxv)
37 
38 !max IF (IFMODEL) CALL SETSHL
39 
40  IF (ifmvbd) THEN
41  ifgeom = .true.
42  IF ( ifflow .AND. .NOT. ifnav ) ifwcno = .true.
43  IF ( ifmelt .AND. .NOT. ifflow ) ifwcno = .true.
44  ENDIF
45 
46  IF (ifflow) THEN
47 
48  ! k call check_cyclic ! fow now; set in .rea file
49 
50  ifield = 1
51  DO iel=1,nelv
52  DO ifc=1,nface
53  cb = cbc(ifc,iel,ifield)
54  CALL chknord(ifalgn,ifnorx,ifnory,ifnorz,ifc,iel)
55  IF ( .NOT. ifstrs ) CALL chkcbc(cb,iel,ifc,ifalgn)
56  IF (cb == 'O ' .OR. cb == 'o ' .OR. &
57  cb == 'ON ' .OR. cb == 'on ' .OR. &
58  cb == 'S ' .OR. cb == 's ' .OR. &
59  cb == 'SL ' .OR. cb == 'sl ' .OR. &
60  cb == 'MM ' .OR. cb == 'mm ' .OR. &
61  cb == 'MS ' .OR. cb == 'ms ') THEN
62  ifvcor = .false.
63  ifeppm(ifc,iel) = .true.
64  ENDIF
65  IF (cb == 'VL ' .OR. cb == 'vl ' .OR. &
66  cb == 'WSL' .OR. cb == 'wsl' .OR. &
67  cb == 'SL ' .OR. cb == 'sl ' .OR. &
68  cb == 'SHL' .OR. cb == 'shl' .OR. &
69  cb == 'MM ' .OR. cb == 'mm ' .OR. &
70  cb == 'MS ' .OR. cb == 'ms ' .OR. &
71  cb == 'O ' .OR. cb == 'o ' .OR. &
72  cb == 'ON ' .OR. cb == 'on ') THEN
73  ifqinp(ifc,iel) = .true.
74  ENDIF
75  IF (cb == 'MS ' .OR. cb == 'ms ' .OR. &
76  cb == 'MM ' .OR. cb == 'mm ' .OR. &
77  cb == 'MSI' .OR. cb == 'msi' ) THEN
78  ifsurt = .true.
79  ENDIF
80  IF (cb == 'WS ' .OR. cb == 'ws ' .OR. &
81  cb == 'WSL' .OR. cb == 'wsl') THEN
82  ifswall = .true.
83  ifcwuz = .true.
84  ENDIF
85  END DO
86  enddo
87  ENDIF
88 
89  IF (ifheat) THEN
90 
91  DO ifield=2,nfield
92  DO iel=1,nelfld(ifield)
93  DO ifc=1,nface
94  cb=cbc(ifc,iel,ifield)
95  IF (cb == 'r ' .OR. cb == 'R ') THEN
96  ifnonl(ifield) = .true.
97  ENDIF
98  enddo
99  enddo
100  END DO
101 
102  ENDIF
103 
104 !max if (ifmhd) call set_ifbcor
105 
106  IF (nhis > 0) THEN
107  iq = 0
108  DO ih=1, nhis
109  IF ( hcode(10,ih) == 'I' ) THEN
110  ifintq = .true.
111  iobj = lochis(1,ih)
112  iq = iq + 1
113  IF (iobj > nobj .OR. iobj < 0) THEN
114  WRITE (6,*) &
115  'ERROR : Undefined Object for integral',iq
116  call exitt
117  ENDIF
118  ENDIF
119  END DO
120  ENDIF
121 
122 ! Establish global consistency of LOGICALS amongst all processors.
123 
124  CALL gllog(ifvcor , .false. )
125  CALL gllog(ifsurt , .true. )
126  CALL gllog(ifswall, .true. )
127  CALL gllog(ifcwuz , .true. )
128  CALL gllog(ifwcno , .true. )
129  DO ifield=2,nfield
130  CALL gllog(ifnonl(ifield), .true. )
131  END DO
132 
133  IF (nid == 0) THEN
134  WRITE (6,*) 'IFTRAN =',iftran
135  WRITE (6,*) 'IFFLOW =',ifflow
136  WRITE (6,*) 'IFHEAT =',ifheat
137  WRITE (6,*) 'IFSPLIT =',ifsplit
138  WRITE (6,*) 'IFLOMACH =',iflomach
139  WRITE (6,*) 'IFUSERVP =',ifuservp
140  WRITE (6,*) 'IFUSERMV =',ifusermv
141  WRITE (6,*) 'IFSTRS =',ifstrs
142  WRITE (6,*) 'IFCHAR =',ifchar
143  WRITE (6,*) 'IFCYCLIC =',ifcyclic
144  WRITE (6,*) 'IFAXIS =',ifaxis
145  WRITE (6,*) 'IFMVBD =',ifmvbd
146  WRITE (6,*) 'IFMELT =',ifmelt
147  WRITE (6,*) 'IFMODEL =',ifmodel
148  WRITE (6,*) 'IFKEPS =',ifkeps
149  WRITE (6,*) 'IFMOAB =',ifmoab
150  WRITE (6,*) 'IFNEKNEK =',ifneknek
151  WRITE (6,*) 'IFSYNC =',ifsync
152  WRITE (6,*) ' '
153  WRITE (6,*) 'IFVCOR =',ifvcor
154  WRITE (6,*) 'IFINTQ =',ifintq
155  WRITE (6,*) 'IFCWUZ =',ifcwuz
156  WRITE (6,*) 'IFSWALL =',ifswall
157  WRITE (6,*) 'IFGEOM =',ifgeom
158  WRITE (6,*) 'IFSURT =',ifsurt
159  WRITE (6,*) 'IFWCNO =',ifwcno
160  DO ifield=1,nfield
161  WRITE (6,*) ' '
162  WRITE (6,*) 'IFTMSH for field',ifield,' = ',iftmsh(ifield)
163  WRITE (6,*) 'IFADVC for field',ifield,' = ',ifadvc(ifield)
164  WRITE (6,*) 'IFNONL for field',ifield,' = ',ifnonl(ifield)
165  END DO
166  WRITE (6,*) ' '
167  if (param(99) > 0) write(6,*) 'Dealiasing enabled, lxd=', lxd
168  ENDIF
169 
170  RETURN
171 END SUBROUTINE setlog
172 
173 !-----------------------------------------------------------------------
176 SUBROUTINE chknord (IFALGN,IFNORX,IFNORY,IFNORZ,IFC,IEL)
177  use kinds, only : dp
178  use size_m, only : ndim, nx1, ny1
179  use geom, only : unx, uny, unz
180  implicit none
181 
182  LOGICAL :: IFALGN,IFNORX,IFNORY,IFNORZ
183  integer :: ifc, iel
184 
185  integer :: ix, iy, ncpf
186  real(DP) :: sumx, sumy, sumz, tolnor
187 
188  sumx = 0.0
189  sumy = 0.0
190  sumz = 0.0
191  tolnor = 1.0e-3
192  ifalgn = .false.
193  ifnorx = .false.
194  ifnory = .false.
195  ifnorz = .false.
196 
197  IF (ndim == 2) THEN
198 
199  ncpf = nx1
200  DO ix=1,nx1
201  sumx = sumx + abs( abs(unx(ix,1,ifc,iel)) - 1.0 )
202  sumy = sumy + abs( abs(uny(ix,1,ifc,iel)) - 1.0 )
203  END DO
204  sumx = sumx / ncpf
205  sumy = sumy / ncpf
206  IF ( sumx < tolnor ) THEN
207  ifnorx = .true.
208  ifalgn = .true.
209  ENDIF
210  IF ( sumy < tolnor ) THEN
211  ifnory = .true.
212  ifalgn = .true.
213  ENDIF
214 
215  ELSE
216 
217  ncpf = nx1*nx1
218  DO ix=1,nx1
219  DO iy=1,ny1
220  sumx = sumx + abs( abs(unx(ix,iy,ifc,iel)) - 1.0 )
221  sumy = sumy + abs( abs(uny(ix,iy,ifc,iel)) - 1.0 )
222  sumz = sumz + abs( abs(unz(ix,iy,ifc,iel)) - 1.0 )
223  enddo
224  END DO
225  sumx = sumx / ncpf
226  sumy = sumy / ncpf
227  sumz = sumz / ncpf
228  IF ( sumx < tolnor ) THEN
229  ifnorx = .true.
230  ifalgn = .true.
231  ENDIF
232  IF ( sumy < tolnor ) THEN
233  ifnory = .true.
234  ifalgn = .true.
235  ENDIF
236  IF ( sumz < tolnor ) THEN
237  ifnorz = .true.
238  ifalgn = .true.
239  ENDIF
240 
241  ENDIF
242 
243  RETURN
244 END SUBROUTINE chknord
245 
246 !-----------------------------------------------------------------------
247 SUBROUTINE chkaxcb()
248  use size_m, only : ndim, nelv
249  use input, only : cbc
250  implicit none
251 
252  CHARACTER(3) CB
253  integer :: ifld, nface, iel, ifc
254 
255  ifld = 1
256  nface = 2*ndim
257 
258  DO iel=1,nelv
259  DO ifc=1,nface
260  cb = cbc(ifc,iel,ifld)
261  IF (cb == 'A ' .AND. ifc /= 1) goto 9000
262  enddo
263  END DO
264 
265  RETURN
266 
267  9000 WRITE (6,*) ' Element face on the axis of symmetry must be FACE 1'
268  WRITE (6,*) ' Element',iel,' face',ifc,' is on the axis.'
269  call exitt
270 
271 END SUBROUTINE chkaxcb
272 !-----------------------------------------------------------------------
274 SUBROUTINE chkcbc (CB,IEL,IFC,IFALGN)
275  implicit none
276  CHARACTER(3) CB
277  LOGICAL :: IFALGN
278  integer :: iel, ifc
279 
280 ! Laplacian formulation only
281 
282  IF (cb == 'SH ' .OR. cb == 'sh ' .OR. &
283  cb == 'SHL' .OR. cb == 'shl' .OR. &
284  cb == 'S ' .OR. cb == 's ' .OR. &
285  cb == 'SL ' .OR. cb == 'sl ' .OR. &
286  cb == 'MM ' .OR. cb == 'mm ' .OR. &
287  cb == 'MS ' .OR. cb == 'ms ' .OR. &
288  cb == 'MSI' .OR. cb == 'msi' ) goto 9001
289  IF ( .NOT. ifalgn .AND. &
290  (cb == 'ON ' .OR. cb == 'on ' .OR. cb == 'SYM') ) goto 9010
291  RETURN
292 
293  9001 WRITE (6,*) ' Illegal traction boundary conditions detected for'
294  goto 9999
295  9010 WRITE (6,*) ' Mixed B.C. on a side nonaligned with either the X,Y, &
296  &or Z axis detected for'
297  9999 WRITE (6,*) ' Element',iel,' side',ifc,'.'
298  WRITE (6,*) ' Selected option only allowed for STRESS FORMULATION'
299  WRITE (6,*) ' Execution terminates'
300  call exitt
301 END SUBROUTINE chkcbc
302 !-----------------------------------------------------------------------
304 SUBROUTINE bcmask
305  use kinds, only : dp
306  use size_m, only : nx1, ny1, nz1, nelv, ndim, nfield
307  use input, only : ifmvbd, ifflow, cbc, ifstrs, ifheat, ipscal, ifmhd
308  use input, only : ifaxis
309  use soln, only : v1mask, v2mask, v3mask, pmask, omask, tmask
310  use tstep, only : ifield, nelfld
311  implicit none
312 
313  character(3) :: cb
314  character(1) :: cb1(3)
315  equivalence(cb1,cb)
316 
317  logical :: ifalgn,ifnorx,ifnory,ifnorz
318 
319  integer :: nfaces, nxyz, nel, ntot, iel, iface
320 
321  nfaces=2*ndim
322  nxyz =nx1*ny1*nz1
323 
324 
325 ! Masks for moving mesh
326 
327  IF (ifmvbd) THEN
328 #if 0
329  ifield = 0
330  CALL stsmask(w1mask,w2mask,w3mask)
331  do e=1,nelv
332  do f=1,nfaces
333  if (cbc(f,e,1) == 'msi' .OR. cbc(f,e,1) == 'msi') then
334  call facev(w1mask,e,f,0.0,nx1,ny1,nz1)
335  call facev(w2mask,e,f,0.0,nx1,ny1,nz1)
336  call facev(w3mask,e,f,0.0,nx1,ny1,nz1)
337  endif
338  enddo
339  enddo
340 #endif
341  ENDIF
342 
343 ! Masks for flow variables
344 
345  IF (ifflow) THEN
346  ifield = 1
347  nel = nelfld(ifield)
348  ntot = nxyz*nel
349 
350  ! Pressure mask
351 
352  pmask => tmask(:,:,:,:,1)
353 #if 0
354  DO iel=1,nelv
355  DO iface=1,nfaces
356  cb=cbc(iface,iel,ifield)
357  IF (cb == 'O ' .OR. cb == 'ON ') &
358  CALL facev(pmask,iel,iface,0.0,nx1,ny1,nz1)
359  enddo
360  END DO
361 
362  ! Zero out mask at Neumann-Dirichlet interfaces
363 
364  CALL dsop(pmask,'MUL')
365 #endif
366 
367  ! Velocity masks
368 
369  ! write(6,*) 'MASK ifstrs',ifstrs,ifield
370  ! call exitt
371  IF (ifstrs) THEN
372 !max CALL STSMASK (V1MASK,V2MASK,V3MASK)
373  ELSE
374 
375  v1mask = 1._dp
376 #if 1
377  v2mask = 1._dp
378  v3mask = 1._dp
379 #else
380  v2mask => v1mask
381  v3mask => v1mask
382 #endif
383 !max if (ifaxis) CALL RONE( OMASK,NTOT)
384 
385  DO iel=1,nelv
386  DO iface=1,nfaces
387  cb =cbc(iface,iel,ifield)
388  CALL chknord(ifalgn,ifnorx,ifnory,ifnorz,iface,iel)
389 
390  ! All-Dirichlet boundary conditions
391 
392  IF (cb == 'v ' .OR. cb == 'V ' .OR. cb == 'vl ' .OR. &
393  cb == 'VL ' .OR. cb == 'W ') THEN
394  CALL facev(v1mask,iel,iface,0.0,nx1,ny1,nz1)
395  CALL facev(v2mask,iel,iface,0.0,nx1,ny1,nz1)
396  CALL facev(v3mask,iel,iface,0.0,nx1,ny1,nz1)
397  cycle
398  ENDIF
399 
400  ! Mixed-Dirichlet-Neumann boundary conditions
401 
402  IF (cb == 'SYM') THEN
403  IF ( .NOT. ifalgn .OR. ifnorx ) &
404  CALL facev(v1mask,iel,iface,0.0,nx1,ny1,nz1)
405  IF ( ifnory ) &
406  CALL facev(v2mask,iel,iface,0.0,nx1,ny1,nz1)
407  IF ( ifnorz ) &
408  CALL facev(v3mask,iel,iface,0.0,nx1,ny1,nz1)
409  cycle
410  ENDIF
411  IF (cb == 'ON ') THEN
412  write(*,*) "Oops! vmasks aren't the same"
413 #if 0
414  IF ( ifnory .OR. ifnorz ) &
415  CALL facev(v1mask,iel,iface,0.0,nx1,ny1,nz1)
416  IF ( .NOT. ifalgn .OR. ifnorx .OR. ifnorz ) &
417  CALL facev(v2mask,iel,iface,0.0,nx1,ny1,nz1)
418  IF ( .NOT. ifalgn .OR. ifnorx .OR. ifnory ) &
419  CALL facev(v3mask,iel,iface,0.0,nx1,ny1,nz1)
420  cycle
421 #endif
422  ENDIF
423  IF (cb == 'A ') THEN
424  write(*,*) "Oops! vmasks aren't the same"
425 #if 0
426  CALL facev(v2mask,iel,iface,0.0,nx1,ny1,nz1)
427  CALL facev( omask,iel,iface,0.0,nx1,ny1,nz1)
428 #endif
429  ENDIF
430  enddo
431  END DO
432 
433  if (ifaxis) CALL dsop( omask,'MUL')
434  call opdsop(v1mask,v2mask,v3mask,'MUL') ! no rotation for mul
435 
436 
437  ENDIF
438 
439  ENDIF
440 
441 ! Masks for passive scalars +
442 ! k and e if k-e turbulence modem:
443 ! k = nfield-1
444 ! e = nfield
445 
446  IF (ifheat) THEN
447 
448  DO ifield=2,nfield
449  ipscal = ifield-1
450  nel = nelfld(ifield)
451  ntot = nxyz*nel
452  tmask(:,:,:,:,ipscal) = 1._dp
453 #if 0
454  DO iel=1,nel
455  DO iface=1,nfaces
456  cb =cbc(iface,iel,ifield)
457 
458  ! Assign mask values.
459 
460  IF (cb == 'T ' .OR. cb == 't ' .OR. &
461  (cb == 'A ' .AND. ifaziv) .OR. &
462  cb == 'MCI' .OR. cb == 'MLI' .OR. &
463  cb == 'KD ' .OR. cb == 'kd ' .OR. &
464  cb == 'ED ' .OR. cb == 'ed ' .OR. &
465  cb == 'KW ' .OR. cb == 'KWS' .OR. cb == 'EWS') &
466  CALL facev(tmask(1,1,1,1,ipscal), &
467  iel,iface,0.0,nx1,ny1,nz1)
468  enddo
469  END DO
470  CALL dsop(tmask(1,1,1,1,ipscal),'MUL')
471 #endif
472  END DO
473 
474  ENDIF
475 
476 ! Masks for B-field
477 
478  if (ifmhd) then
479  write(*,*) "Oops: ifmhd"
480 #if 0
481  ifield = ifldmhd
482  nel = nelfld(ifield)
483  ntot = nxyz*nel
484 
485  ! B-field pressure mask
486 
487  call rone(bpmask,ntot)
488  do iel=1,nelv
489  do iface=1,nfaces
490  cb=cbc(iface,iel,ifield)
491  if (cb == 'O ' .OR. cb == 'ON ') &
492  call facev(bpmask,iel,iface,0.0,nx1,ny1,nz1)
493  enddo
494  enddo
495 
496  ! Zero out mask at Neumann-Dirichlet interfaces
497 
498  call dsop(bpmask,'MUL')
499 
500  ! B-field masks
501 
502  if (ifstrs) then
503 !max call stsmask (b1mask,b2mask,b3mask)
504  else
505 
506  call rone(b1mask,ntot)
507  call rone(b2mask,ntot)
508  call rone(b3mask,ntot)
509 
510  do iel=1,nelv
511  do iface=1,nfaces
512  cb =cbc(iface,iel,ifield)
513  call chknord(ifalgn,ifnorx,ifnory,ifnorz,iface,iel)
514 
515  if (cb == 'v ' .OR. cb == 'V ' .OR. cb == 'vl ' .OR. &
516  cb == 'VL ' .OR. cb == 'W ') then
517 
518  ! All-Dirichlet boundary conditions
519 
520  call facev(b1mask,iel,iface,0.0,nx1,ny1,nz1)
521  call facev(b2mask,iel,iface,0.0,nx1,ny1,nz1)
522  call facev(b3mask,iel,iface,0.0,nx1,ny1,nz1)
523 
524  elseif (cb == 'SYM') then
525 
526  ! Mixed-Dirichlet-Neumann boundary conditions
527 
528  if ( .NOT. ifalgn .OR. ifnorx ) &
529  call facev(b1mask,iel,iface,0.0,nx1,ny1,nz1)
530  if ( ifnory ) &
531  call facev(b2mask,iel,iface,0.0,nx1,ny1,nz1)
532  if ( ifnorz ) &
533  call facev(b3mask,iel,iface,0.0,nx1,ny1,nz1)
534 
535  elseif (cb == 'ON ') then
536 
537  ! Mixed-Dirichlet-Neumann boundary conditions
538 
539  if ( ifnory .OR. ifnorz ) &
540  call facev(b1mask,iel,iface,0.0,nx1,ny1,nz1)
541  if ( .NOT. ifalgn .OR. ifnorx .OR. ifnorz ) &
542  call facev(b2mask,iel,iface,0.0,nx1,ny1,nz1)
543  if ( .NOT. ifalgn .OR. ifnorx .OR. ifnory ) &
544  call facev(b3mask,iel,iface,0.0,nx1,ny1,nz1)
545 
546  elseif (cb == 'A ') then
547 
548  ! axisymmetric centerline
549 
550  call facev(b2mask,iel,iface,0.0,nx1,ny1,nz1)
551 
552  else
553 
554  if ( cb1(1) == 'd' ) &
555  call facev(b1mask,iel,iface,0.0,nx1,ny1,nz1)
556  if ( cb1(2) == 'd' ) &
557  call facev(b2mask,iel,iface,0.0,nx1,ny1,nz1)
558  if ( cb1(3) == 'd' .AND. if3d ) &
559  call facev(b3mask,iel,iface,0.0,nx1,ny1,nz1)
560 
561  endif
562  enddo
563  enddo
564 
565  call dsop(b1mask,'MUL')
566  call dsop(b2mask,'MUL')
567  if (ndim == 3) call dsop(b3mask,'MUL')
568  endif
569 #endif
570  endif
571 
572  RETURN
573 END SUBROUTINE bcmask
574 !-----------------------------------------------------------------------
588 SUBROUTINE bcdirvc(V1,V2,V3,mask1,mask2,mask3)
589  use kinds, only : dp
590  use size_m, only : nx1, ny1, nz1, ndim, lelv
591  use ctimer, only : icalld, tusbc, nusbc, etime1, dnekclock
592  use input, only : if3d, cbc, bc, ifstrs
593  use tstep, only : ifield, nelfld
594  implicit none
595 
596  REAL(DP) :: V1(nx1,ny1,nz1,lelv),V2(nx1,ny1,nz1,lelv) &
597  ,V3(nx1,ny1,nz1,lelv)
598  real(DP) :: mask1(nx1,ny1,nz1,lelv),mask2(nx1,ny1,nz1,lelv) &
599  ,mask3(nx1,ny1,nz1,lelv)
600 
601  character(3) cb
602  character(1) :: cb1(3)
603  equivalence(cb1,cb)
604 
605  logical :: ifonbc
606 
607  integer :: nfaces, nxyz, nel, ntot, isweep, ie, iface
608  real(DP) :: bc1, bc2, bc3
609 
610  ifonbc = .false.
611 
612 #ifndef NOTIMER
613  if (icalld == 0) then
614  tusbc=0.0
615  nusbc=0
616  icalld=icalld+1
617  endif
618  nusbc=nusbc+1
619  etime1=dnekclock()
620 #endif
621 
622 
623  nfaces=2*ndim
624  nxyz =nx1*ny1*nz1
625  nel =nelfld(ifield)
626  ntot =nxyz*nel
627 
628 #if 0
629  allocate(tmp1(nx1,ny1,nz1,nelfld(ifield)) &
630  , tmp2(nx1,ny1,nz1,nelfld(ifield)) &
631  , tmp3(nx1,ny1,nz1,nelfld(ifield)) )
632 
633  tmp1 = 0._dp; tmp2 = 0._dp
634  IF (if3d) tmp3 = 0._dp
635 #endif
636 
637 ! Velocity boundary conditions
638 
639 ! write(6,*) 'BCDIRV: ifield',ifield
640  DO isweep=1,2
641  DO ie=1,nel
642  DO iface=1,nfaces
643  cb = cbc(iface,ie,ifield)
644  bc1 = bc(1,iface,ie,ifield)
645  bc2 = bc(2,iface,ie,ifield)
646  bc3 = bc(3,iface,ie,ifield)
647 
648  IF (cb == 'V ' .OR. cb == 'VL ' .OR. &
649  cb == 'WS ' .OR. cb == 'WSL') THEN
650  write(*,*) "Oops: CB = v, vl, ws, wsl"
651 #if 0
652  CALL facev(tmp1,ie,iface,bc1,nx1,ny1,nz1)
653  CALL facev(tmp2,ie,iface,bc2,nx1,ny1,nz1)
654  IF (if3d) CALL facev(tmp3,ie,iface,bc3,nx1,ny1,nz1)
655  IF ( ifqinp(iface,ie) ) &
656  CALL globrot(tmp1(1,1,1,ie),tmp2(1,1,1,ie), &
657  tmp3(1,1,1,ie),ie,iface)
658 #endif
659  ENDIF
660 
661  IF (cb == 'v ' .OR. cb == 'vl ' .OR. &
662  cb == 'ws ' .OR. cb == 'wsl' .OR. &
663  cb == 'mv ' .OR. cb == 'mvn' .OR. &
664  cb1(1) == 'd' .OR. cb1(2) == 'd' .OR. cb1(3) == 'd') then
665  write(*,*) "Oops: CB = something bad"
666 #if 0
667  call faceiv(cb,tmp1(1,1,1,ie),tmp2(1,1,1,ie), &
668  tmp3(1,1,1,ie),ie,iface,nx1,ny1,nz1)
669 
670  IF ( ifqinp(iface,ie) ) &
671  CALL globrot(tmp1(1,1,1,ie),tmp2(1,1,1,ie), &
672  tmp3(1,1,1,ie),ie,iface)
673 #endif
674  ENDIF
675 
676  IF (cb == 'ON ' .OR. cb == 'on ') then ! 5/21/01 pff
677  write(*,*) "Oops: CB = ON"
678 #if 0
679  ifonbc = .true.
680  CALL faceiv('v ',tmp1(1,1,1,ie),tmp2(1,1,1,ie), &
681  tmp3(1,1,1,ie),ie,iface,nx1,ny1,nz1)
682 #endif
683  ENDIF
684  enddo
685  END DO
686 #if 0
687  DO ie=1,nel
688  DO iface=1,nfaces
689  IF (cbc(iface,ie,ifield) == 'W ') THEN
690  CALL facev(tmp1,ie,iface,0.0,nx1,ny1,nz1)
691  CALL facev(tmp2,ie,iface,0.0,nx1,ny1,nz1)
692  IF (if3d) CALL facev(tmp3,ie,iface,0.0,nx1,ny1,nz1)
693  ENDIF
694  enddo
695  END DO
696 
697  ! Take care of Neumann-Dirichlet shared edges...
698 
699  if (isweep == 1) then
700  call opdsop(tmp1,tmp2,tmp3,'MXA')
701  else
702  call opdsop(tmp1,tmp2,tmp3,'MNA')
703  endif
704 #endif
705  END DO
706 ! Copy temporary array to velocity arrays.
707 
708  IF ( .NOT. ifstrs ) THEN
709  v1 = v1 * mask1
710  v2 = v2 * mask2
711  IF (if3d) v3 = v3 * mask3
712  if (ifonbc) then
713  write(*,*) "Oops: ifonbc"
714 #if 0
715  call antimsk1(tmp1,mask1,ntot)
716  call antimsk1(tmp2,mask2,ntot)
717  if (if3d) call antimsk1(tmp3,mask3,ntot)
718 #endif
719  endif
720  ELSE
721  write(*,*) "Oops: ifstrs"
722 #if 0
723  IF (ifmodel) THEN
724  CALL copy(tmq1,tmp1,ntot)
725  CALL copy(tmq2,tmp2,ntot)
726  IF (ndim == 3) CALL copy(tmq3,tmp3,ntot)
727  CALL amask(tmp1,tmp2,tmp3,tmq1,tmq2,tmq3,nelv)
728  ENDIF
729  CALL rmask(v1,v2,v3,nelv)
730 #endif
731  ENDIF
732 
733 #if 0
734  v1 = v1 + tmp1; v2 = v2 + tmp2
735  IF (if3d) v3 = v3 + tmp3
736 #endif
737 
738 
739 #ifndef NOTIMER
740  tusbc=tusbc+(dnekclock()-etime1)
741 #endif
742 
743  RETURN
744 END SUBROUTINE bcdirvc
745 !-----------------------------------------------------------------------
748 SUBROUTINE bcdirsc(S)
749  use kinds, only : dp
750  use size_m, only : lx1, ly1, lz1, lelt
751  use size_m, only : nx1, ny1, nz1, ndim, nfield
752  use ctimer, only : icalld, tusbc, nusbc, etime1, dnekclock
753  use input, only : cbc, bc
754  use soln, only : tmask
755  use tstep, only : ifield, nelfld
756  implicit none
757 
758  real(DP) :: S(lx1,ly1,lz1,lelt)
759  real(DP), allocatable :: tmp(:,:,:,:)
760 
761  CHARACTER(3) CB
762 
763  integer :: ifld, nfaces, nxyz, nel, ntot, isweep, ie, iface, nfldt
764  real(DP) :: BC1, BC2, BC3, BC4, BCK, BCE
765 
766 #ifndef NOTIMER
767  if (icalld == 0) then
768  tusbc=0.0
769  nusbc=0
770  icalld=icalld+1
771  endif
772  nusbc=nusbc+1
773  etime1=dnekclock()
774 #endif
775 
776  ifld = 1
777  nfaces = 2*ndim
778  nxyz = nx1*ny1*nz1
779  nel = nelfld(ifield)
780  ntot = nxyz*nel
781  nfldt = nfield - 1
782 
783  allocate(tmp(nx1,ny1,nz1,nelfld(ifield)))
784  tmp = 0._dp
785 
786 ! Temperature boundary condition
787 
788  DO isweep=1,2
789 
790 #if 0
791  IF (ifmodel .AND. ifkeps .AND. ifield >= nfldt) &
792  CALL turbwbc(tmp,tma,smu)
793 #endif
794 
795  DO ie=1,nel
796  DO iface=1,nfaces
797  cb=cbc(iface,ie,ifield)
798  bc1=bc(1,iface,ie,ifield)
799  bc2=bc(2,iface,ie,ifield)
800  bc3=bc(3,iface,ie,ifield)
801  bc4=bc(4,iface,ie,ifield)
802  bck=bc(4,iface,ie,ifld)
803  bce=bc(5,iface,ie,ifld)
804  IF (cb == 'T ') CALL facev(tmp,ie,iface,bc1,nx1,ny1,nz1)
805  IF (cb == 'MCI') CALL facev(tmp,ie,iface,bc4,nx1,ny1,nz1)
806  IF (cb == 'MLI') CALL facev(tmp,ie,iface,bc4,nx1,ny1,nz1)
807  IF (cb == 'KD ') CALL facev(tmp,ie,iface,bck,nx1,ny1,nz1)
808  IF (cb == 'ED ') CALL facev(tmp,ie,iface,bce,nx1,ny1,nz1)
809  IF (cb == 't ' .OR. cb == 'kd ' .OR. cb == 'ed ') then
810  write(*,*) "Oops: CB = t, kd, or ed"
811 !max CALL FACEIS (CB,TMP(1,1,1,IE),IE,IFACE,NX1,NY1,NZ1)
812  ENDIF
813  END DO
814  enddo
815 
816  ! Take care of Neumann-Dirichlet shared edges...
817 
818  IF (isweep == 1) CALL dsop(tmp,'MXA')
819  IF (isweep == 2) CALL dsop(tmp,'MNA')
820  END DO
821 
822 ! Copy temporary array to temperature array.
823 
824  s = s * tmask(:,:,:,:,ifield-1) + tmp
825 
826 #ifndef NOTIMER
827  tusbc=tusbc+(dnekclock()-etime1)
828 #endif
829 
830  RETURN
831 END SUBROUTINE bcdirsc
832 
833 !-----------------------------------------------------------------------
840 SUBROUTINE bcneusc(S,ITYPE)
841  use kinds, only : dp
842  use size_m, only : lx1, ly1, lz1, lelt, nx1, ny1, nz1, ndim
843  use ctimer, only : icalld, tusbc, nusbc, etime1, dnekclock
844  use nekuse, only : hc, tinf, hrad, flux
845  use geom, only : area
846  use input, only : cbc, bc
847  use geom, only : bm1
848  use parallel, only : lglel
849  use soln, only : t
850  use tstep, only : ifield, nelfld
851  implicit none
852 
853  real(DP), intent(out) :: S(lx1,ly1,lz1,lelt)
854  integer :: itype
855 
856  CHARACTER(3) CB
857 
858  real(DP) :: ts
859  integer :: nfaces, nxyz, nel, ntot
860  integer :: ie, iface, ieg, ia, ix, iy, iz
861  integer :: kx1, kx2, ky1, ky2, kz1, kz2
862 
863 #ifndef NOTIMER
864  if (icalld == 0) then
865  tusbc=0.0
866  nusbc=0
867  icalld=icalld+1
868  endif
869  nusbc=nusbc+1
870  etime1=dnekclock()
871 #endif
872 
873  nfaces=2*ndim
874  nxyz =nx1*ny1*nz1
875  nel =nelfld(ifield)
876  ntot =nxyz*nel
877  s = 0._dp
878 
879  IF (itype == -1) THEN
880 
881  ! Compute diagonal contributions to accomodate Robin boundary conditions
882 
883  DO ie=1,nel
884  DO iface=1,nfaces
885  ieg=lglel(ie)
886  cb =cbc(iface,ie,ifield)
887  IF (cb == 'C ' .OR. cb == 'c ' .OR. &
888  cb == 'R ' .OR. cb == 'r ') THEN
889 
890  IF (cb == 'C ') hc = bc(2,iface,ie,ifield)
891  IF (cb == 'R ') THEN
892  tinf = bc(1,iface,ie,ifield)
893  hrad = bc(2,iface,ie,ifield)
894  ENDIF
895  ia=0
896 
897  ! IA is areal counter, assumes advancing fastest index first. (IX...IY...IZ)
898 
899  CALL facind(kx1,kx2,ky1,ky2,kz1,kz2,nx1,ny1,nz1,iface)
900  DO iz=kz1,kz2
901  DO iy=ky1,ky2
902  DO ix=kx1,kx2
903  ia = ia + 1
904  ts = t(ix,iy,iz,ie,ifield-1)
905  IF (cb == 'c ' .OR. cb == 'r ') THEN
906  CALL nekasgn(ix,iy,iz,ie)
907  CALL userbc(ix,iy,iz,iface,ieg)
908  ENDIF
909  IF (cb == 'r ' .OR. cb == 'R ') &
910  hc = hrad * (tinf**2 + ts**2) * (tinf + ts)
911  s(ix,iy,iz,ie) = s(ix,iy,iz,ie) + &
912  hc*area(ia,1,iface,ie)/bm1(ix,iy,iz,ie)
913  enddo
914  enddo
915  END DO
916  ENDIF
917  enddo
918  END DO
919  ENDIF
920  IF (itype == 1) THEN
921 
922  ! Add passive scalar fluxes to rhs
923 
924  DO ie=1,nel
925  DO iface=1,nfaces
926  ieg=lglel(ie)
927  cb =cbc(iface,ie,ifield)
928  IF (cb == 'F ' .OR. cb == 'f ' .OR. &
929  cb == 'C ' .OR. cb == 'c ' .OR. &
930  cb == 'R ' .OR. cb == 'r ' ) THEN
931 
932  IF (cb == 'F ') flux=bc(1,iface,ie,ifield)
933  IF (cb == 'C ') flux=bc(1,iface,ie,ifield) &
934  *bc(2,iface,ie,ifield)
935  IF (cb == 'R ') THEN
936  tinf=bc(1,iface,ie,ifield)
937  hrad=bc(2,iface,ie,ifield)
938  ENDIF
939 
940  ! Add local weighted flux values to rhs, S.
941 
942  ! IA is areal counter, assumes advancing fastest index first. (IX...IY...IZ)
943  ia=0
944  CALL facind(kx1,kx2,ky1,ky2,kz1,kz2,nx1,ny1,nz1,iface)
945  DO iz=kz1,kz2
946  DO iy=ky1,ky2
947  DO ix=kx1,kx2
948  ia = ia + 1
949  ts = t(ix,iy,iz,ie,ifield-1)
950  IF (cb == 'f ') THEN
951  CALL nekasgn(ix,iy,iz,ie)
952  CALL userbc(ix,iy,iz,iface,ieg)
953  ENDIF
954  IF (cb == 'c ') THEN
955  CALL nekasgn(ix,iy,iz,ie)
956  CALL userbc(ix,iy,iz,iface,ieg)
957  flux = tinf*hc
958  ENDIF
959  IF (cb == 'r ') THEN
960  CALL nekasgn(ix,iy,iz,ie)
961  CALL userbc(ix,iy,iz,iface,ieg)
962  ENDIF
963  IF (cb == 'R ' .OR. cb == 'r ') &
964  flux = hrad*(tinf**2 + ts**2)*(tinf + ts) * tinf
965 
966  ! Add computed fluxes to boundary surfaces:
967 
968  s(ix,iy,iz,ie) = s(ix,iy,iz,ie) &
969  + flux*area(ia,1,iface,ie)
970  enddo
971  enddo
972  enddo
973  ENDIF
974  enddo
975  END DO
976  ENDIF
977 
978 #ifndef NOTIMER
979  tusbc=tusbc+(dnekclock()-etime1)
980 #endif
981 
982  RETURN
983 END SUBROUTINE bcneusc
984 !-----------------------------------------------------------------------
1014 SUBROUTINE nekasgn (IX,IY,IZ,IEL)
1015  use size_m
1016  use geom
1017  use input
1018  use nekuse
1019  use soln
1020  use tstep
1021  implicit none
1022 
1023  integer, intent(in) :: ix, iy, iz, iel
1024 
1025  CHARACTER(3) CB
1026 
1027  integer :: ips
1028 
1029  x = xm1(ix,iy,iz,iel)
1030  y = ym1(ix,iy,iz,iel)
1031  z = zm1(ix,iy,iz,iel)
1032  r = x**2+y**2
1033  IF (r > 0.0) r=sqrt(r)
1034  IF (x /= 0.0 .OR. y /= 0.0) theta = atan2(y,x)
1035 
1036  ux = vx(ix,iy,iz,iel)
1037  uy = vy(ix,iy,iz,iel)
1038  uz = vz(ix,iy,iz,iel)
1039  temp = t(ix,iy,iz,iel,1)
1040  DO 100 ips=1,npscal
1041  ps(ips) = t(ix,iy,iz,iel,ips+1)
1042  100 END DO
1043  udiff = vdiff(ix,iy,iz,iel,ifield)
1044  utrans= vtrans(ix,iy,iz,iel,ifield)
1045 
1046  cbu = cb
1047 
1048  RETURN
1049 END SUBROUTINE nekasgn
1050 !-----------------------------------------------------------------------
1051 SUBROUTINE lfalse (IFA,N)
1052  implicit none
1053  LOGICAL :: IFA(*)
1054  integer :: n, i
1055  DO 100 i=1,n
1056  ifa(i)= .false.
1057  100 END DO
1058  RETURN
1059 END SUBROUTINE lfalse
1060 !-----------------------------------------------------------------------
1061 SUBROUTINE unitvec (X,Y,Z,N)
1062  use kinds, only : dp
1063  implicit none
1064  integer, intent(in) :: n
1065  real(DP), intent(inout) :: X(n),Y(n),Z(n)
1066  real(DP) :: xlngth
1067  integer :: i
1068  DO i=1,n
1069  xlngth = sqrt( x(i)**2 + y(i)**2 + z(i)**2 )
1070  IF (xlngth /= 0.0) THEN
1071  x(i) = x(i)/xlngth
1072  y(i) = y(i)/xlngth
1073  z(i) = z(i)/xlngth
1074  ENDIF
1075  END DO
1076  RETURN
1077 END SUBROUTINE unitvec
1078 !-----------------------------------------------------------------------
tstep
cleaned
Definition: tstep_mod.F90:2
chknord
subroutine chknord(IFALGN, IFNORX, IFNORY, IFNORZ, IFC, IEL)
Check direction of normal of an element face for alignment with the X, Y, or Z axis.
Definition: bdry.F90:176
input
Input parameters from preprocessors.
Definition: input_mod.F90:11
bcneusc
subroutine bcneusc(S, ITYPE)
Apply Neumann boundary conditions to surface of scalar, S. Use IFIELD as a guide to which boundary co...
Definition: bdry.F90:840
soln
Definition: soln_mod.F90:1
bcdirvc
subroutine bcdirvc(V1, V2, V3, mask1, mask2, mask3)
Apply Dirichlet boundary conditions to surface of vector (V1,V2,V3). Use IFIELD as a guide to which b...
Definition: bdry.F90:588
exitt
void exitt()
Definition: comm_mpi.F90:411
setlog
subroutine setlog()
Subroutine to initialize logical flags.
Definition: bdry.F90:3
chkaxcb
subroutine chkaxcb()
Definition: bdry.F90:247
ctimer::dnekclock
real(dp) function dnekclock()
Definition: ctimer_mod.F90:103
parallel::lglel
integer function lglel(iel)
Definition: parallel_mod.F90:176
copy
subroutine copy(a, b, n)
Definition: math.F90:52
bcdirsc
subroutine bcdirsc(S)
Apply Dirichlet boundary conditions to surface of scalar, S. Use IFIELD as a guide to which boundary ...
Definition: bdry.F90:748
chkcbc
subroutine chkcbc(CB, IEL, IFC, IFALGN)
Check for illegal boundary conditions.
Definition: bdry.F90:274
dsop
subroutine dsop(u, op)
generalization of dssum to other reducers.
Definition: dssum.F90:134
parallel
cleaned
Definition: parallel_mod.F90:2
turbo
cleaned
Definition: turbo_mod.F90:2
lfalse
subroutine lfalse(IFA, N)
Definition: bdry.F90:1051
geom
Geometry arrays.
Definition: geom_mod.F90:2
nekasgn
subroutine nekasgn(IX, IY, IZ, IEL)
Assign NEKTON variables for definition (by user) of boundary conditions at collocation point (IX...
Definition: bdry.F90:1014
bcmask
subroutine bcmask
Zero out masks corresponding to Dirichlet boundary points.
Definition: bdry.F90:304
unitvec
subroutine unitvec(X, Y, Z, N)
Definition: bdry.F90:1061
nekuse
cleaned
Definition: nekuse_mod.F90:2
y
static double y[NR *NS *NT *N]
Definition: obbox_test.c:31
facev
subroutine facev(a, ie, iface, val, nx, ny, nz)
Assign the value VAL to face(IFACE,IE) of array A. IFACE is the input in the pre-processor ordering s...
Definition: connect1.F90:1041
gllog
subroutine gllog(la, lb)
If ANY LA=LB, then ALL LA=LB.
Definition: math.F90:402
facind
subroutine facind(kx1, kx2, ky1, ky2, kz1, kz2, nx, ny, nz, iface)
ifcase in preprocessor notation
Definition: connect1.F90:1019
z
static double z[NR *NS *NT *N]
Definition: obbox_test.c:31
Generated on Tue Jul 14 2015 10:52:58 for Nek5000 by   doxygen 1.8.8