-
Notifications
You must be signed in to change notification settings - Fork 4
/
utilities-2P.usr
731 lines (658 loc) · 21.6 KB
/
utilities-2P.usr
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
c-----------------------------------------------------------------------
subroutine area_average(phi,loc,coord,norm,phia)
C
C Compute area average of phi() on the
C plane normal to norm() with intercept coord = loc
C
include 'SIZE'
include 'TOTAL'
integer ielem,iside,i,i0,i1,j,j0,j1,k,k0,k1
real phi(1),wrt(1),loc,coord(1),norm(3),phia
real fnorm(3),dp,a1,phia1
loc=get_nearest_face(loc,coord,norm)
phia=0.0
phia1=0.0
do ielem=1,nelv
do iside=1,ndim*2
call facind(i0,i1,j0,j1,k0,k1,lx1,ly1,lz1,iside)
i=(i0+i1)/2 !just use the point in the middle of the face
j=(j0+j1)/2
k=(k0+k1)/2
ipoint=i+(j-1)*lx1+(k-1)*lx1*ly1+(ielem-1)*lx1*ly1*lz1
if(abs(coord(ipoint)-loc).lt.1.0d-8)then
call getSnormal(fnorm,i,j,k,iside,ielem)
dp=fnorm(1)*norm(1)+fnorm(2)*norm(2)
if(if3d) dp=dp+fnorm(3)*norm(3)
if(abs(1.0d0-dp).lt.1.0d-8) then
do i=i0,i1
do j=j0,j1
do k=k0,k1
ipoint=i+(j-1)*lx1+(k-1)*lx1*ly1+(ielem-1)*lx1*ly1*lz1
if ((iside.eq.1).or.(iside.eq.3)) then
a1=area(i,k,iside,ielem)
elseif((iside.eq.2).or.(iside.eq.4)) then
a1=area(j,k,iside,ielem)
else
a1=area(i,j,iside,ielem)
endif
phia=phia+phi(ipoint)*a1
phia1=phia1+a1
enddo
enddo
enddo
endif
endif
enddo
enddo
phia=glsum(phia,1)
phia1=glsum(phia1,1)
phia=phia/phia1
return
end
C-----------------------------------------------------------------------
subroutine count_boundaries
include 'SIZE'
include 'TOTAL'
integer lxyz,ielem,iside,n
parameter(lxyz=lx1*ly1*lz1)
character*3 uid(ldimt+1)
integer wcount(ldimt+1),symcount(ldimt+1),ocount(ldimt+1)
integer tcount(ldimt+1),fcount(ldimt+1),axicount(ldimt+1)
integer inscount(ldimt+1),pcount(ldimt+1),othcount(ldimt+1)
integer vcount(ldimt+1),trcount(ldimt+1),ukncount(ldimt+1)
integer vreacount(ldimt+1),treacount(ldimt+1)
call izero(wcount,ldimt+1)
call izero(trcount,ldimt+1)
call izero(vcount,ldimt+1)
call izero(vreacount,ldimt+1)
call izero(symcount,ldimt+1)
call izero(ocount,ldimt+1)
call izero(tcount,ldimt+1)
call izero(treacount,ldimt+1)
call izero(axicount,ldimt+1)
call izero(inscount,ldimt+1)
call izero(fcount,ldimt+1)
call izero(pcount,ldimt+1)
call izero(othcount,ldimt+1)
call izero(ukncount,ldimt+1)
do ifld=1,nfield
do ielem=1,nelv
do iside=1,2*ldim
if(cbc(iside,ielem,ifld).eq.'W ')then
wcount(ifld)=wcount(ifld)+1
elseif(cbc(iside,ielem,ifld).eq.'shl')then
trcount(ifld)=trcount(ifld)+1
elseif(cbc(iside,ielem,ifld).eq.'v ')then
vcount(ifld)=vcount(ifld)+1
elseif(cbc(iside,ielem,ifld).eq.'V ')then
vreacount(ifld)=vreacount(ifld)+1
elseif(cbc(iside,ielem,ifld).eq.'t ')then
tcount(ifld)=tcount(ifld)+1
elseif(cbc(iside,ielem,ifld).eq.'T ')then
treacount(ifld)=treacount(ifld)+1
elseif(cbc(iside,ielem,ifld).eq.'O ')then
ocount(ifld)=ocount(ifld)+1
elseif(cbc(iside,ielem,ifld).eq.'P ')then
pcount(ifld)=pcount(ifld)+1
elseif(cbc(iside,ielem,ifld).eq.'f ')then
fcount(ifld)=fcount(ifld)+1
elseif(cbc(iside,ielem,ifld).eq.'I ')then
inscount(ifld)=inscount(ifld)+1
elseif(cbc(iside,ielem,ifld).eq.'SYM')then
symcount(ifld)=symcount(ifld)+1
elseif(cbc(iside,ielem,ifld).eq.'A ')then
axicount(ifld)=axicount(ifld)+1
elseif(cbc(iside,ielem,ifld).ne.'E ')then
if(ukncount(ifld).eq.0) then !handle one unknown BC
uid(ifld)=cbc(iside,ielem,ifld)
ukncount(ifld)=1
elseif(cbc(iside,ielem,ifld).eq.uid(ifld)) then
ukncount(ifld)=ukncount(ifld)+1
else
othcount(ifld)=othcount(ifld)+1 !multiple unknown BCs
endif
endif
enddo
enddo
wcount(ifld)=iglsum(wcount(ifld),1)
trcount(ifld)=iglsum(trcount(ifld),1)
vcount(ifld)=iglsum(vcount(ifld),1)
vreacount(ifld)=iglsum(vreacount(ifld),1)
tcount(ifld)=iglsum(tcount(ifld),1)
treacount(ifld)=iglsum(treacount(ifld),1)
ocount(ifld)=iglsum(ocount(ifld),1)
pcount(ifld)=iglsum(pcount(ifld),1)
fcount(ifld)=iglsum(fcount(ifld),1)
inscount(ifld)=iglsum(inscount(ifld),1)
symcount(ifld)=iglsum(symcount(ifld),1)
axicount(ifld)=iglsum(axicount(ifld),1)
othcount(ifld)=iglsum(othcount(ifld),1)
ukncount(ifld)=iglsum(ukncount(ifld),1)
enddo
if(nid.eq.0) then
write(*,*)
write(*,255) 'Found the following Boundary Conditions'
write(*,*)
do ifld=1,nfield
write(*,254) 'for field',ifld,':'
if(wcount(ifld).gt.0)write(*,256)'Wall',wcount(ifld)
if(trcount(ifld).gt.0)write(*,256)'Traction',trcount(ifld)
if(vcount(ifld).gt.0)write(*,256)'Velocity',vcount(ifld)
if(vreacount(ifld).gt.0)write(*,256)'Velocity (REA)',
& vreacount(ifld)
if(tcount(ifld).gt.0)write(*,256)'Dirichlet',tcount(ifld)
if(treacount(ifld).gt.0)write(*,256)'Dirichlet (REA)',
& treacount(ifld)
if(pcount(ifld).gt.0)write(*,256)'Periodic',pcount(ifld)
if(fcount(ifld).gt.0)write(*,256)'Flux',fcount(ifld)
if(ocount(ifld).gt.0)write(*,256)'Outlet',ocount(ifld)
if(inscount(ifld).gt.0)write(*,256)'Insulated',inscount(ifld)
if(symcount(ifld).gt.0)write(*,256)'Symmetry',symcount(ifld)
if(axicount(ifld).gt.0)write(*,256)
& 'Axisymmetric',axicount(ifld)
if(ukncount(ifld).gt.0)write(*,257)uid(ifld),ukncount(ifld)
if(othcount(ifld).gt.0)write(*,256)'Other',othcount(ifld)
write(*,*)
enddo
endif
254 format(5x,a,i2,a)
255 format(2x,a)
256 format(2x,a16,1x,i6)
257 format(2x,'Unknown boundary of type: ',a,1x,i6)
return
end
C-----------------------------------------------------------------------
subroutine div_check(phi)
real phi(1)
if(phi(1).ne.phi(1)) call exitt
return
end
C-----------------------------------------------------------------------
subroutine get_limits(phi,phimin,phimax,phiave,dphi,phip,rmsphi,n)
include 'SIZE'
include 'TOTAL'
integer lxyze
parameter(lxyze=lx1*ly1*lz1*lelv)
real phi(1),phip(1),phimin,phimax,phiave,dphi,rmsphi
integer i,n,ntot
rmsphi=0.0
dphi=0.0
if(istep.ge.1) then
do i=1,n
dphi=max(dphi,abs(phip(i)-phi(i)))
rmsphi=rmsphi+bm1(i,1,1,1)*(phip(i)-phi(i))**2
enddo
rmsphi=glsum(rmsphi,1)/volvm1
rmsphi=sqrt(rmsphi)/dt
endif
dphi=glmax(dphi,1)
dphi=dphi/dt
phimin=glmin(phi,n)
phimax=glmax(phi,n)
phiave=glsc2(phi,bm1,n)/volvm1
return
end
c-----------------------------------------------------------------------
subroutine get_limits_nodt(phi,phimin,phimax,phiave,n)
implicit none
include 'SIZE'
include 'TOTAL'
integer n
real phi(1),phimin,phimax,phiave
real glmin,glmax,glsc2
phimin=glmin(phi,n)
phimax=glmax(phi,n)
phiave=glsc2(phi,bm1,n)/volvm1
return
end
c-----------------------------------------------------------------------
real function get_nearest(loc,coord)
include 'mpif.h'
include 'SIZE'
include 'TOTAL'
integer ipoint,ierr
real loc,coord(1),ds(2),dsg(2)
ds(1)=1.0d30
do ipoint=1,lx1*ly1*lz1*nelv
if(abs(coord(ipoint)-loc).lt.ds(1)) then
ds(1)=abs(coord(ipoint)-loc)
ds(2)=coord(ipoint)
endif
enddo
call MPI_ALLREDUCE(ds,dsg,1,MPI_2DOUBLE_PRECISION,MPI_MINLOC
& ,MPI_COMM_WORLD,ierr)
get_nearest=dsg(2)
return
end
c-----------------------------------------------------------------------
real function get_nearest_face(loc,coord,norm)
include 'mpif.h'
include 'SIZE'
include 'TOTAL'
integer ielem,iside,i0,i1,j0,j1,k0,k1,i,j,k,ierr
real loc,coord(1),ds(2),dsg(2),norm(3),fnorm(3),dp
ds(1)=1.0d30
do ielem=1,nelv
do iside=1,ldim*2
call facind(i0,i1,j0,j1,k0,k1,lx1,ly1,lz1,iside)
i=(i0+i1)/2
j=(j0+j1)/2
k=(k0+k1)/2
ipoint=i+(j-1)*lx1+(k-1)*lx1*ly1+lx1*ly1*lz1*(ielem-1)
call getSnormal(fnorm,i,j,k,iside,ielem)
dp=fnorm(1)*norm(1)+fnorm(2)*norm(2)
if(if3d) dp=dp+fnorm(3)*norm(3)
if((1.0d0-abs(dp)).lt.1.0d-8)then
if(abs(coord(ipoint)-loc).lt.ds(1)) then
ds(1)=abs(coord(ipoint)-loc)
ds(2)=coord(ipoint)
endif
endif
enddo
enddo
call MPI_ALLREDUCE(ds,dsg,1,MPI_2DOUBLE_PRECISION,MPI_MINLOC
& ,MPI_COMM_WORLD,ierr)
get_nearest_face=dsg(2)
return
end
C-----------------------------------------------------------------------
subroutine max_y_p(wd)
include 'SIZE'
include 'TOTAL'
real wd(nx1,ny1,nz1,nelv)
integer ielem,iside,i,i0,i1,j,j0,j1,k,k0,k1
integer ne,mask(lx1,ly1,lz1)
real dwmax,dwmin
ne=nx1*ny1*nz1
dwmax=-1.0d30
dwmin=1.0d30
do ielem=1,nelv
call izero(mask,ne) !mask all points in the element
do iside=1,2*ndim !unmask all points just off the wall
if(cbc(iside,ielem,1).eq.'W ') then
call backpts(i0,i1,j0,j1,k0,k1,iside)
do i=i0,i1
do j=j0,j1
do k=k0,k1
mask(i,j,k)=1
enddo
enddo
enddo
endif
enddo
do iside=1,2*ndim !remask all points on the wall
if(cbc(iside,ielem,1).eq.'W ') then
call facind(i0,i1,j0,j1,k0,k1,nx1,ny1,nz1,iside)
do k=k0,k1
do j=j0,j1
do i=i0,i1
mask(i,j,k)=0
enddo
enddo
enddo
endif
enddo
do i=1,ne
if(mask(i,1,1).ne.0)then
dwmax=max(dwmax,wd(i,1,1,ielem))
dwmin=min(dwmin,wd(i,1,1,ielem))
endif
enddo
enddo
dwmax=glmax(dwmax,1)
dwmin=glmin(dwmin,1)
if(nid.eq.0) then
write(*,256) 'maximum y_p = ',dwmax
write(*,256) 'minimum y_p = ',dwmin
write(*,*)
endif
256 format(2x,a,es15.5)
return
end
C-----------------------------------------------------------------------
subroutine weighted_average(phi,wrt,loc,coord,norm,phia)
C
C Compute planar averages of phi() weighted by wrt() on the
C plane normal to norm() with intercept coord = loc
C
include 'SIZE'
include 'TOTAL'
integer ielem,iside,i,i0,i1,j,j0,j1,k,k0,k1
real phi(1),wrt(1),loc,coord(1),norm(3),phia
real fnorm(3),dp,a1,phia1
loc=get_nearest_face(loc,coord,norm)
phia=0.0
phia1=0.0
do ielem=1,nelv
do iside=1,ndim*2
call facind(i0,i1,j0,j1,k0,k1,lx1,ly1,lz1,iside)
i=(i0+i1)/2 !just use the point in the middle of the face
j=(j0+j1)/2
k=(k0+k1)/2
ipoint=i+(j-1)*lx1+(k-1)*lx1*ly1+(ielem-1)*lx1*ly1*lz1
if(abs(coord(ipoint)-loc).lt.1.0d-8)then
call getSnormal(fnorm,i,j,k,iside,ielem)
dp=fnorm(1)*norm(1)+fnorm(2)*norm(2)
if(if3d) dp=dp+fnorm(3)*norm(3)
if(abs(1.0d0-dp).lt.1.0d-8) then
do i=i0,i1
do j=j0,j1
do k=k0,k1
ipoint=i+(j-1)*lx1+(k-1)*lx1*ly1+(ielem-1)*lx1*ly1*lz1
if ((iside.eq.1).or.(iside.eq.3)) then
a1=area(i,k,iside,ielem)
elseif((iside.eq.2).or.(iside.eq.4)) then
a1=area(j,k,iside,ielem)
else
a1=area(i,j,iside,ielem)
endif
phia=phia+phi(ipoint)*wrt(ipoint)*a1
phia1=phia1+wrt(ipoint)*a1
enddo
enddo
enddo
endif
endif
enddo
enddo
phia=glsum(phia,1)
phia1=glsum(phia1,1)
phia=phia/phia1
return
end
C-----------------------------------------------------------------------
subroutine planar_average_weighted(phia,phi,wrt,w1,w2)
include 'SIZE'
include 'GEOM'
include 'PARALLEL'
include 'WZ'
include 'ZPER'
real phi(nx1*ny1,nz1,nelv),wrt(nx1*ny1,nz1,nelv),phia(nz1,nelz)
real w1(nz1,nelz),w2(nz1,nelz) !work arrays
integer e,eg,ez,melxy,nz,i,k
real zz,aa
melxy=nelx*nely !number of elements in the plane
nz=nz1*nelz !number of z-slices
if(melxy.lt.1)then
if(nio.eq.0)write(*,256)'nelx*nely'
return
elseif(nelz.lt.1) then
if(nio.eq.0)write(*,256)'nelz'
return
elseif(melxy.gt.lelx*lely) then
if(nio.eq.0)write(*,257)'nelx*nely','lelx*lely'
return
elseif(nelz.gt.lelz) then
if(nio.eq.0)write(*,257)'nelz','lelz'
return
endif
256 format(5x,'ERROR: ',a,' must be at least 1!')
257 format(5x,'ERROR: ',a,' must be less than ',a,'!')
call rzero(phia,nz)
call rzero(w1,nz)
do e=1,nelt
eg=lglel(e)
ez=1+(eg-1)/melxy !z-slice id
do k=1,nz1
do i=1,nx1*ny1
zz=(1.0-zgm1(k,3))/2.0
aa=zz*area(i,1,5,e)+(1.0-zz)*area(i,1,6,e)
w1(k,ez)=w1(k,ez)+aa*wrt(i,k,e)
phia(k,ez)=phia(k,ez)+aa*wrt(i,k,e)*phi(i,k,e)
enddo
if(abs(w1(k,ez)).lt.1.0d-10) w1(k,ez)=1.0d-10
enddo
enddo
call gop(phia,w2,'+ ',nz)
call gop(w1,w2,'+ ',nz)
call invcol2(phia,w1,nz)
return
end
C-----------------------------------------------------------------------
subroutine x_planar_average(phia,phi,w1,w2)
include 'SIZE'
include 'GEOM'
include 'PARALLEL'
include 'WZ'
include 'ZPER'
real phi(nx1,ny1,nz1,nelv),wrt(nx1,ny1,nz1,nelv),phia(nx1,nelx)
real w1(nx1,nelx),w2(nx1,nelx) !work arrays
integer e,eg,ex,nx,i,j,estride
real xx,aa
if(ldim.gt.2) then
write(*,'(5x,a)')
& "x-average routine only written for 2D genbox meshes!"
return
endif
nx=nx1*nelx !number of z-slices
call rzero(phia,nx)
call rzero(w1,nx)
do e=1,nelt
eg=lglel(e)
ex=mod(eg,nelx) !x-slice id
if(ex.eq.0)ex=nelx
do i=1,nx1
do j=1,ny1
xx=(1.0-zgm1(i,1))/2.0
aa=zz*area(j,1,4,e)+(1.0-zz)*area(j,1,2,e)
w1(i,ex)=w1(i,ex)+aa
phia(i,ex)=phia(i,ex)+aa*phi(i,j,1,e)
enddo
enddo
enddo
call gop(phia,w2,'+ ',nx)
call gop(w1,w2,'+ ',nx)
call invcol2(phia,w1,nx)
return
end
C-----------------------------------------------------------------------
subroutine x_average_weighted(phia,phi,wrt,w1,w2)
include 'SIZE'
include 'GEOM'
include 'PARALLEL'
include 'WZ'
include 'ZPER'
real phi(nx1,ny1,nz1,nelv),wrt(nx1,ny1,nz1,nelv),phia(nx1,nelx)
real w1(nx1,nelx),w2(nx1,nelx) !work arrays
integer e,eg,ex,nx,i,j,estride
real xx,aa
if(lz1.gt.1) then
write(*,'(5x,a)')
& "x-average routine only written for 2D meshes!"
return
endif
nx=nx1*nelx !number of z-slices
call rzero(phia,nx)
call rzero(w1,nx)
do e=1,nelt
eg=lglel(e)
ex=mod(eg,nelx) !x-slice id
if(ex.eq.0)ex=nelx
do i=1,nx1
do j=1,ny1
xx=(1.0-zgm1(i,1))/2.0
aa=zz*area(j,1,4,e)+(1.0-zz)*area(j,1,2,e)
w1(i,ex)=w1(i,ex)+aa*wrt(i,j,1,e)
phia(i,ex)=phia(i,ex)+aa*wrt(i,j,1,e)*phi(i,j,1,e)
enddo
if(abs(w1(i,ex)).lt.1.0d-10) w1(i,ex)=1.0d-10
enddo
enddo
call gop(phia,w2,'+ ',nx)
call gop(w1,w2,'+ ',nx)
call invcol2(phia,w1,nx)
return
end
C-----------------------------------------------------------------------
subroutine y_p_limits(wd,ypmin,ypmax,ypave,utmin,utmax,utave,iph)
implicit none
include 'SIZE'
include 'TOTAL'
C
C NOTE: min value should work if domain has internal corners
C ...hasn't been tested
integer e,i,i0,i1,j,j0,j1,k,k0,k1,iw,jw,kw,i2,j2
integer ipt,wpt,estrd,isd,jsd,iph
integer imsk(lx1,ly1,lz1,lelv)
real gradu(lx1*ly1*lz1,3,3),wd(1)
real tau(3),norm(3),vsca,tauw,utau,rho,mu
real ypmin,ypmax,yp,ypave,vol,utmin,utmax,utave
real glmin,glmax,glsum
logical flag, idid
data idid /.false./
save idid, imsk
ypmin=1.0d30
ypmax=-1.0d30
ypave=0.0
utmin=1.0d30
utmax=-1.0d30
utave=0.0
vol=0.0
C first build the mask
if(.not.idid)then
idid=.true.
call izero(imsk,nx1*ny1*nz1*nelv)
do e=1,nelv
do isd=1,2*ndim
if(cbc(isd,e,1).eq.'W ') then
call backpts(i0,i1,j0,j1,k0,k1,isd)
do k=k0,k1
do j=j0,j1
do i=i0,i1
imsk(i,j,k,e)=1
enddo
enddo
enddo
endif
enddo
do isd=1,2*ndim
if(cbc(isd,e,1).eq.'W ') then
call facind(i0,i1,j0,j1,k0,k1,lx1,ly1,lz1,isd)
do k=k0,k1
do j=j0,j1
do i=i0,i1
imsk(i,j,k,e)=0
enddo
enddo
enddo
endif
enddo
enddo
endif
do e=1,nelv
flag=.true.
do isd=1,2*ndim
if(cbc(isd,e,1).eq.'W ')then
estrd=(e-1)*nx1*ny1*nz1
if(flag)then
call gradm11(gradu(1,1,1),gradu(1,1,2),gradu(1,1,3)
& ,vkx(1,1,1,1,iph),e)
call gradm11(gradu(1,2,1),gradu(1,2,2),gradu(1,2,3)
& ,vky(1,1,1,1,iph),e)
if(if3d)
& call gradm11(gradu(1,3,1),gradu(1,3,2),gradu(1,3,3)
& ,vkz(1,1,1,1,iph),e)
flag=.false.
endif
call backpts(i0,i1,j0,j1,k0,k1,isd)
do k=k0,k1
do j=j0,j1
do i=i0,i1
if(imsk(i,j,k,e).eq.1) then
iw=i
jw=j
kw=k
if (isd.eq.1) then
jw=1
elseif(isd.eq.2) then
iw=nx1
elseif(isd.eq.3) then
jw=ny1
elseif(isd.eq.4) then
iw=1
elseif(isd.eq.5) then
kw=1
else
kw=nx1
endif
call getSnormal(norm,iw,jw,kw,isd,e)
ipt=i +(j -1)*nx1+(k -1)*nx1*ny1
wpt=iw+(jw-1)*nx1+(kw-1)*nx1*ny1
mu=vdiffk(iw,jw,kw,e,iph)
rho=vtransk(iw,jw,kw,e,iph)
do i2=1,ldim
tau(i2)=0.0
do j2=1,ldim
tau(i2)=tau(i2)+
& mu*(gradu(wpt,i2,j2)+gradu(wpt,j2,i2))*norm(j2)
enddo
enddo
vsca=0.0
do i2=1,ldim
vsca=vsca+tau(i2)*norm(i2)
enddo
tauw=0.0
do i2=1,ldim
tauw=tauw+(tau(i2)-vsca*norm(i2))**2
enddo
tauw=sqrt(tauw)
utau=sqrt(tauw/rho)
yp=wd(ipt+estrd)*utau*rho/mu
ypmin=min(ypmin,yp)
ypmax=max(ypmax,yp)
utave=utave+utau*bm1(i,j,k,e)
utmin=min(utau,utmin)
utmax=max(utau,utmax)
ypave=ypave+yp*bm1(i,j,k,e)
vol=vol+bm1(i,j,k,e)
endif
enddo
enddo
enddo
endif
enddo
enddo
ypmin=glmin(ypmin,1)
ypmax=glmax(ypmax,1)
ypave=glsum(ypave,1)
utmin=glmin(utmin,1)
utmax=glmax(utmax,1)
utave=glsum(utave,1)
vol=glsum(vol,1)
ypave=ypave/vol
utave=utave/vol
return
end
c-----------------------------------------------------------------------
subroutine backpts(i0,i1,j0,j1,k0,k1,isd)
implicit none
include 'SIZE'
integer i0,i1,j0,j1,k0,k1,isd
i0=1
j0=1
k0=1
i1=nx1
j1=ny1
k1=nz1
if(isd.eq.1) then
j0=2
j1=2
elseif(isd.eq.2) then
i0=nx1-1
i1=nx1-1
elseif(isd.eq.3) then
j0=ny1-1
j1=ny1-1
elseif(isd.eq.4) then
i0=2
i1=2
elseif(isd.eq.5) then
k0=2
k1=2
elseif(isd.eq.6) then
k0=nz1-1
k1=nz1-1
endif
return
end
c-----------------------------------------------------------------------