source: CIVL/examples/fortran/nek5000/core/multimesh.f

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100755
File size: 32.7 KB
Line 
1c-----------------------------------------------------------------------
2c Routines for multidomain (neknek) simulations.
3c
4c References:
5c
6c "A spectrally accurate method for overlapping grid solution of
7c incompressible Navier–Stokes equations" Brandon E. Merrill,
8c Yulia T. Peet, Paul F. Fischer, and James W. Lottes, J. Comp. Phys.
9c 307 (2016) 60-93.
10c
11c "Stability analysis of interface temporal discretization in grid
12c overlapping methods," Y. Peet, P.F. Fischer, SIAM J. Numer. Anal.
13c 50 (6) (2012) 3375–3401.
14c-----------------------------------------------------------------------
15 subroutine setup_neknek_wts
16
17 include 'SIZE'
18 include 'TOTAL'
19
20 ntot1 = lx1*ly1*lz1*nelt
21c Initialize unity partition function to 1
22 call rone(upf,ntot1)
23 call col3(bm1ms,bm1,upf,ntot1)
24
25 return
26 end
27c-------------------------------------------------------------
28 subroutine neknek_setup
29
30 include 'SIZE'
31 include 'TOTAL'
32 real dxf,dyf,dzf
33
34 integer icalld
35 save icalld
36 data icalld /0/
37
38 integer nfld_neknek
39 common /inbc/ nfld_neknek
40
41 if (icalld.eq.0.and.nid.eq.0) write(6,*) 'setup neknek'
42
43 if (nsessmax.eq.1)
44 $ call exitti('set nsessmax > 1 in SIZE!$',nsessmax)
45
46 call setup_neknek_wts
47
48 if (icalld.eq.0) then
49 ! just in case we call setup from usrdat2
50 call fix_geom
51 call geom_reset(1)
52
53 call set_intflag
54 call neknekmv
55 if (nid.eq.0) write(6,*) 'session id:', idsess
56 if (nid.eq.0) write(6,*) 'extrapolation order:', ninter
57 if (nid.eq.0) write(6,*) 'nfld_neknek:', nfld_neknek
58
59 nfld_min = iglmin_ms(nfld_neknek,1)
60 nfld_max = iglmax_ms(nfld_neknek,1)
61 if (nfld_min .ne. nfld_max) then
62 nfld_neknek = nfld_min
63 if (nid.eq.0) write(6,*)
64 $ 'WARNING: reset nfld_neknek to ', nfld_neknek
65 endif
66 endif
67
68c Figure out the displacement for the first mesh
69 call setup_int_neknek(dxf,dyf,dzf) !sets up interpolation for 2 meshes
70
71c exchange_points finds the processor and element number at
72c comm_world level and displaces the 1st mesh back
73 call exchange_points(dxf,dyf,dzf)
74
75 if (icalld.eq.0) then
76 if(nio.eq.0) write(6,'(A,/)') ' done :: setup neknek'
77 icalld = 1
78 endif
79
80 return
81 end
82c-------------------------------------------------------------
83 subroutine set_intflag
84 include 'SIZE'
85 include 'TOTAL'
86 include 'NEKNEK'
87 character*3 cb
88 character*2 cb2
89 equivalence (cb2,cb)
90 integer j,e,f
91
92c Set interpolation flag: points with boundary condition = 'int'
93c get intflag=1.
94c
95c Boundary conditions are changed back to 'v' or 't'.
96
97 nfaces = 2*ldim
98
99 nflag=nelt*nfaces
100 call izero(intflag,nflag)
101
102 do j=1,nfield
103 nel = nelfld(j)
104 do e=1,nel
105 do f=1,nfaces
106 cb=cbc(f,e,j)
107 if (cb2.eq.'in') then
108 intflag(f,e)=1
109 if (j.ge.2) cbc(f,e,j)='t '
110 if (j.eq.1) cbc(f,e,j)='v '
111c if (cb.eq.'inp') cbc(f,e,ifield)='on ' ! Pressure
112c if (cb.eq.'inp') cbc(f,e,ifield)='o ' ! Pressure
113 if (cb.eq.'inp') cbc(f,e,j)='o ' ! Pressure
114 endif
115 enddo
116 enddo
117 enddo
118
119c zero out valint
120 do i=1,nfld_neknek
121 call rzero(valint(1,1,1,1,i),lx1*ly1*lz1*nelt)
122 enddo
123
124 return
125 end
126c------------------------------------------------------------------------
127 subroutine bcopy
128 include 'SIZE'
129 include 'TOTAL'
130 include 'NEKNEK'
131 integer k,i,n
132
133 if (.not.ifneknekc) return
134
135 n = lx1*ly1*lz1*nelt
136
137 do k=1,nfld_neknek
138 call copy(bdrylg(1,k,2),bdrylg(1,k,1),n)
139 call copy(bdrylg(1,k,1),bdrylg(1,k,0),n)
140 call copy(bdrylg(1,k,0),valint(1,1,1,1,k),n)
141 enddo
142
143c Order of extrpolation is contolled by the parameter NINTER contained
144c in NEKNEK. First order interface extrapolation, NINTER=1 (time lagging)
145c is activated. It is unconditionally stable. If you want to use
146c higher-order interface extrapolation schemes, you need to increase
147c ngeom to ngeom=3-5 for scheme to be stable.
148
149
150 if (NINTER.eq.1.or.istep.eq.1) then
151 c0=1.
152 c1=0.
153 c2=0.
154 else if (NINTER.eq.2.or.istep.eq.2) then
155 c0=2.
156 c1=-1.
157 c2=0.
158 else
159 c0=3.
160 c1=-3.
161 c2=1.
162 endif
163
164 do k=1,nfld_neknek
165 do i=1,n
166 valint(i,1,1,1,k) =
167 $ c0*bdrylg(i,k,0)+c1*bdrylg(i,k,1)+c2*bdrylg(i,k,2)
168 enddo
169 enddo
170
171 return
172 end
173c---------------------------------------------------------------------
174 subroutine chk_outflow ! Assign neighbor velocity to outflow if
175 ! characteristics are going the wrong way.
176
177 include 'SIZE'
178 include 'TOTAL'
179 include 'NEKUSE'
180 include 'NEKNEK'
181 integer e,eg,f
182
183 nface = 2*ldim
184 do e=1,nelv
185 do f=1,nface
186 if (cbc(f,e,1).eq.'o ') then
187 eg = lglel(e)
188 call facind (i0,i1,j0,j1,k0,k1,lx1,ly1,lz1,f)
189 l=0
190 do k=k0,k1
191 do j=j0,j1
192 do i=i0,i1
193 l=l+1
194 vo=vx(i,j,k,e)*unx(l,1,f,e)
195 $ +vy(i,j,k,e)*uny(l,1,f,e)
196 $ +vz(i,j,k,e)*unz(l,1,f,e)
197 if (vo.lt.0) then ! We have inflow
198 cbu = cbc(f,e,1)
199 call userbc(i,j,k,f,eg)
200 vx(i,j,k,e) = ux
201 vy(i,j,k,e) = uy
202 vz(i,j,k,e) = uz
203 endif
204 enddo
205 enddo
206 enddo
207 endif
208 enddo
209 enddo
210
211 return
212 end
213c-----------------------------------------------------------------------
214 subroutine neknekmv()
215 include 'SIZE'
216 include 'TOTAL'
217
218 integer imove
219
220 imove=1
221 if (ifmvbd) imove=0
222
223 iglmove = iglmin_ms(imove,1)
224
225 if (iglmove.eq.0) then
226 ifneknekm=.true.
227 endif
228
229 return
230 end
231c-----------------------------------------------------------------------
232 subroutine setup_int_neknek(dxf,dyf,dzf)
233 include 'SIZE'
234 include 'TOTAL'
235 include 'NEKUSE'
236 include 'NEKNEK'
237 include 'mpif.h'
238
239 real dx1,dy1,dz1,dxf,dyf,dzf,mx_glob,mn_glob
240 integer i,j,k,n,ntot2,npall
241 common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
242c THIS ROUTINE DISPLACES THE FIRST MESH AND SETUPS THE FINDPTS
243c THE MESH IS DISPLACED BACK TO ORIGINAL POSITION IN EXCH_POINTS
244
245c Get total number of processors and number of p
246 npall = 0
247 do i=1,nsessions
248 npall = npall+npsess(i-1)
249 enddo
250
251c Get diamter of the domain
252 mx_glob = glmax_ms(xm1,lx1*ly1*lz1*nelt)
253 mn_glob = glmin_ms(xm1,lx1*ly1*lz1*nelt)
254 dx1 = mx_glob-mn_glob
255
256 dxf = 10.+dx1
257 dyf = 0.
258 dzf = 0.
259
260c Displace MESH 1
261 ntot = lx1*ly1*lz1*nelt
262 if (idsess.eq.0) then
263 call cadd(xm1,-dxf,ntot)
264 endif
265
266c Setup findpts
267 tol = 5e-13
268 npt_max = 128
269 nxf = 2*lx1 ! fine mesh for bb-test
270 nyf = 2*ly1
271 nzf = 2*lz1
272 bb_t = 0.01 ! relative size to expand bounding boxes by
273
274 if (istep.gt.1) call fgslib_findpts_free(inth_multi2)
275 call fgslib_findpts_setup(inth_multi2,mpi_comm_world,npall,ldim,
276 & xm1,ym1,zm1,lx1,ly1,lz1,
277 & nelt,nxf,nyf,nzf,bb_t,ntot,ntot,
278 & npt_max,tol)
279
280 return
281 end
282c-----------------------------------------------------------------------
283 subroutine exchange_points(dxf,dyf,dzf)
284 include 'SIZE'
285 include 'TOTAL'
286 include 'NEKUSE'
287 include 'NEKNEK'
288 integer i,j,k,n
289 common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
290 integer jsend(nmaxl_nn)
291 common /exchr/ rsend(ldim*nmaxl_nn)
292 integer rcode_all(nmaxl_nn),elid_all(nmaxl_nn)
293 integer proc_all(nmaxl_nn)
294 real dist_all(nmaxl_nn)
295 real rst_all(nmaxl_nn*ldim)
296 integer e,ip,iface,nel,nfaces,ix,iy,iz
297 integer kx1,kx2,ky1,ky2,kz1,kz2,idx,nxyz,nxy
298 integer icalld
299 save icalld
300 data icalld /0/
301
302c Look for boundary points with Diriclet b.c. (candidates for
303c interpolation)
304
305 ifield = 1
306 if (ifheat) ifield = 2
307
308 nfaces = 2*ldim
309 nel = nelfld(ifield)
310
311 nxyz = lx1*ly1*lz1
312 ntot = nxyz*nel
313 nxy = lx1*ly1
314 call izero(imask,ntot)
315
316c Setup arrays of x,y,zs to send to findpts and indices of boundary
317c points in jsend
318 ip = 0
319 if (idsess.eq.0) then
320 dxf = -dxf
321 endif
322 do e=1,nel
323 do iface=1,nfaces
324 if (intflag(iface,e).eq.1) then
325 call facind (kx1,kx2,ky1,ky2,kz1,kz2,lx1,ly1,lz1,iface)
326 do iz=kz1,kz2
327 do iy=ky1,ky2
328 do ix=kx1,kx2
329 call nekasgn (ix,iy,iz,e)
330 ip=ip+1
331 idx = (e-1)*nxyz+(iz-1)*nxy+(iy-1)*lx1+ix
332 jsend(ip) = idx
333 if (if3d) then
334 rsend(ldim*(ip-1)+1)=x-dxf
335 rsend(ldim*(ip-1)+2)=y
336 rsend(ldim*(ip-1)+3)=z
337 else
338 rsend(ldim*(ip-1)+1)=x-dxf
339 rsend(ldim*(ip-1)+2)=y
340 endif
341
342 if (ip.gt.nmaxl_nn) then
343 write(6,*) nid,
344 & ' ABORT: nbp (current ip) too large',ip,nmaxl_nn
345 call exitt
346 endif
347
348 imask(idx,1,1,1)=1
349
350 enddo
351 enddo
352 enddo
353 endif
354 enddo
355 enddo
356 nbp = ip
357
358 call neknekgsync()
359
360c JL's routine to find which points these procs are on
361 call fgslib_findpts(inth_multi2,rcode_all,1,
362 & proc_all,1,
363 & elid_all,1,
364 & rst_all,ldim,
365 & dist_all,1,
366 & rsend(1),ldim,
367 & rsend(2),ldim,
368 & rsend(3),ldim,nbp)
369
370 call neknekgsync()
371
372c Move mesh 1 back to its original position
373 if (idsess.eq.0) then
374 call cadd(xm1,-dxf,lx1*ly1*lz1*nelt)
375 endif
376
377 ip=0
378 icount=0
379 ierror=0
380
381c Make sure rcode_all is fine
382 do 200 i=1,nbp
383
384 if (rcode_all(i).lt.2) then
385
386 if (rcode_all(i).eq.1.and.dist_all(i).gt.1e-02) then
387 if (ldim.eq.2) write(6,*)
388 & 'WARNING: point on boundary or outside the mesh xy[z]d^2: '
389 if (ldim.eq.3) write(6,*)
390 & 'WARNING: point on boundary or outside the mesh xy[z]d^2: '
391 goto 200
392 endif
393 ip=ip+1
394 rcode(ip) = rcode_all(i)
395 elid(ip) = elid_all(i)
396 proc(ip) = proc_all(i)
397 do j=1,ldim
398 rst(ldim*(ip-1)+j) = rst_all(ldim*(i-1)+j)
399 enddo
400 iList(1,ip) = jsend(i)
401
402 endif ! rcode_all
403
404 200 continue
405
406 ipg = iglsum(ip,1)
407 nbpg = iglsum(nbp,1)
408 if (nid.eq.0) write(6,*) ipg,nbpg,'interface points'
409 npoints_nn = ip
410
411 ierror = iglmax_ms(ierror,1)
412 if (ierror.eq.1) call exitt
413
414 return
415 end
416c-----------------------------------------------------------------------
417 subroutine neknek_exchange
418 include 'SIZE'
419 include 'TOTAL'
420 include 'NEKNEK'
421 include 'CTIMER'
422
423 parameter (lt=lx1*ly1*lz1*lelt,lxyz=lx1*ly1*lz1)
424 common /scrcg/ pm1(lt),wk1(lxyz),wk2(lxyz)
425
426 real fieldout(nmaxl_nn,nfldmax_nn)
427 real field(lx1*ly1*lz1*lelt)
428 integer nv,nt,i,j,k,n,ie,ix,iy,iz,idx,ifld
429
430 if (nio.eq.0) write(6,98)
431 $ ' Multidomain data exchange ... ', nfld_neknek
432 98 format(12x,a,i3)
433
434 etime0 = dnekclock_sync()
435 call neknekgsync()
436 etime1 = dnekclock()
437
438 call mappr(pm1,pr,wk1,wk2) ! Map pressure to pm1
439 nv = lx1*ly1*lz1*nelv
440 nt = lx1*ly1*lz1*nelt
441
442c Interpolate using findpts_eval
443 call field_eval(fieldout(1,1),1,vx)
444 call field_eval(fieldout(1,2),1,vy)
445 if (ldim.eq.3) call field_eval(fieldout(1,ldim),1,vz)
446 call field_eval(fieldout(1,ldim+1),1,pm1)
447 if (nfld_neknek.gt.ldim+1) then
448 do i=ldim+2,nfld_neknek
449 call field_eval(fieldout(1,i),1,t(1,1,1,1,i-ldim-1))
450 enddo
451 endif
452
453c Now we can transfer this information to valint array from which
454c the information will go to the boundary points
455 do i=1,npoints_nn
456 idx = iList(1,i)
457 do ifld=1,nfld_neknek
458 valint(idx,1,1,1,ifld)=fieldout(i,ifld)
459 enddo
460 enddo
461
462 call nekgsync()
463 etime = dnekclock() - etime1
464 tsync = etime1 - etime0
465
466 if (nio.eq.0) write(6,99) istep,
467 $ ' done :: Multidomain data exchange',
468 $ etime, etime+tsync
469 99 format(i11,a,1p2e13.4)
470
471 return
472 end
473c--------------------------------------------------------------------------
474 subroutine field_eval(fieldout,fieldstride,fieldin)
475 include 'SIZE'
476 include 'TOTAL'
477 include 'NEKNEK'
478 real fieldout(1)
479 real fieldin(1)
480 integer fieldstride
481
482c Used for findpts_eval of various fields
483 call fgslib_findpts_eval(inth_multi2,fieldout,fieldstride,
484 & rcode,1,
485 & proc,1,
486 & elid,1,
487 & rst,ldim,npoints_nn,
488 & fieldin)
489
490 return
491 end
492c--------------------------------------------------------------------------
493 subroutine nekneksanchk
494 include 'SIZE'
495 include 'TOTAL'
496 include 'NEKNEK'
497
498 if (nfld_neknek.gt.nfldmax_nn) then
499 call exitti('Error: nfld_neknek > nfldmax:$',idsess)
500 endif
501
502 if (nfld_neknek.eq.0)
503 $ call exitti('Error: set nfld_neknek in usrdat. Session:$',idsess)
504
505 return
506 end
507C--------------------------------------------------------------------------
508 subroutine neknek_xfer_fld(u,ui)
509 include 'SIZE'
510 include 'TOTAL'
511 include 'NEKNEK'
512 real fieldout(nmaxl_nn,nfldmax_nn)
513 real u(1),ui(1)
514 integer nv,nt
515
516 call field_eval(fieldout(1,1),1,u)
517
518c Now we can transfer this information to valint array from which
519c the information will go to the boundary points
520 do i=1,npoints_nn
521 idx = iList(1,i)
522 ui(idx)=fieldout(i,1)
523 enddo
524 call neknekgsync()
525
526 return
527 end
528C--------------------------------------------------------------------------
529 subroutine fix_surface_flux
530 include 'SIZE'
531 include 'TOTAL'
532 include 'NEKNEK'
533 integer e,f
534 common /ctmp1/ work(lx1*ly1*lz1*lelt)
535 integer itchk
536 common /idumochk/ itchk
537 integer icalld
538 save icalld
539 data icalld /0/
540c assume that this routine is called at the end of bcdirvc
541c where all the boundary condition data has been read in for
542c velocity.
543 if (icalld.eq.0) then
544 itchk = 0
545 do e=1,nelv
546 do f=1,2*ldim
547 if (cbc(f,e,1).eq.'o '.or.cbc(f,e,1).eq.'O ') then
548 if (intflag(f,e).eq.0) then
549 itchk = 1
550 endif
551 endif
552 enddo
553 enddo
554 itchk = iglmax(itchk,1)
555 icalld = 1
556 endif
557 if (itchk.eq.1) return
558 dqg=0
559 aqg=0
560 do e=1,nelv
561 do f=1,2*ldim
562 if (cbc(f,e,1).eq.'v '.or.cbc(f,e,1).eq.'V ') then
563 call surface_flux_area(dq,aq
564 $ ,vx,vy,vz,e,f,work)
565 dqg = dqg+dq
566 if (intflag(f,e).eq.1) aqg = aqg+aq
567 endif
568 enddo
569 enddo
570 dqg=glsum(dqg,1) ! sum over all processors for this session
571 aqg=glsum(aqg,1) ! sum over all processors for this session
572 gamma = 0.
573 if (aqg.gt.0) gamma = -dqg/aqg
574c if (nid.eq.0) write(6,104) idsess,istep,time,dqg,aqg,gamma
575c 104 format(i4,i10,1p4e13.4,' NekNek_bdry_flux')
576 do e=1,nelv
577 do f=1,2*ldim
578 if (intflag(f,e).eq.1) then
579 call facind (i0,i1,j0,j1,k0,k1,lx1,ly1,lz1,f)
580 l=0
581 do k=k0,k1
582 do j=j0,j1
583 do i=i0,i1
584 l=l+1
585 vx(i,j,k,e) = vx(i,j,k,e) + gamma*unx(l,1,f,e)
586 vy(i,j,k,e) = vy(i,j,k,e) + gamma*uny(l,1,f,e)
587 if (ldim.eq.3)
588 $ vz(i,j,k,e) = vz(i,j,k,e) + gamma*unz(l,1,f,e)
589 enddo
590 enddo
591 enddo
592 endif
593 enddo
594 enddo
595 return
596 end
597c-----------------------------------------------------------------------
598 subroutine surface_flux_area(dq,aq,qx,qy,qz,e,f,w)
599 include 'SIZE'
600 include 'GEOM'
601 include 'INPUT'
602 include 'PARALLEL'
603 include 'TOPOL'
604 parameter (l=lx1*ly1*lz1)
605 real qx(l,1),qy(l,1),qz(l,1),w(lx1,ly1,lz1)
606 integer e,f
607 call faccl3 (w,qx(1,e),unx(1,1,f,e),f)
608 call faddcl3 (w,qy(1,e),uny(1,1,f,e),f)
609 if (if3d) call faddcl3 (w,qz(1,e),unz(1,1,f,e),f)
610 call dsset(lx1,ly1,lz1)
611 iface = eface1(f)
612 js1 = skpdat(1,iface)
613 jf1 = skpdat(2,iface)
614 jskip1 = skpdat(3,iface)
615 js2 = skpdat(4,iface)
616 jf2 = skpdat(5,iface)
617 jskip2 = skpdat(6,iface)
618 dq = 0
619 aq = 0
620 i = 0
621 do 100 j2=js2,jf2,jskip2
622 do 100 j1=js1,jf1,jskip1
623 i = i+1
624 dq = dq + area(i,1,f,e)*w(j1,j2,1)
625 aq = aq + area(i,1,f,e)
626 100 continue
627 return
628 end
629c-----------------------------------------------------------------------
630 subroutine rescale_x_ms (x,x0,x1)
631 include 'SIZE'
632 real x(1)
633
634 n = lx1*ly1*lz1*nelt
635 xmin = glmin(x,n)
636 xmax = glmax(x,n)
637 xming = glmin_ms(x,n)
638 xmaxg = glmax_ms(x,n)
639
640 if (xmax.le.xmin) return
641
642 scale = (x1-x0)/(xmaxg-xming)
643 x0n = x0 + scale*(xmin-xming)
644
645
646 do i=1,n
647 x(i) = x0n + scale*(x(i)-xmin)
648 enddo
649
650 return
651 end
652c-----------------------------------------------------------------------
653c-----------------------------------------------------------------------
654 subroutine vol_flow_ms
655c
656c
657c Adust flow volume at end of time step to keep flow rate fixed by
658c adding an appropriate multiple of the linear solution to the Stokes
659c problem arising from a unit forcing in the X-direction. This assumes
660c that the flow rate in the X-direction is to be fixed (as opposed to Y-
661c or Z-) *and* that the periodic boundary conditions in the X-direction
662c occur at the extreme left and right ends of the mesh.
663c
664c pff 6/28/98
665c
666 include 'SIZE'
667 include 'TOTAL'
668c
669c Swap the comments on these two lines if you don't want to fix the
670c flow rate for periodic-in-X (or Z) flow problems.
671c
672 parameter (kx1=lx1,ky1=ly1,kz1=lz1,kx2=lx2,ky2=ly2,kz2=lz2)
673c
674 common /cvflow_a/ vxc(kx1,ky1,kz1,lelv)
675 $ , vyc(kx1,ky1,kz1,lelv)
676 $ , vzc(kx1,ky1,kz1,lelv)
677 $ , prc(kx2,ky2,kz2,lelv)
678 $ , vdc(kx1*ky1*kz1*lelv,2)
679 common /cvflow_r/ flow_rate,base_flow,domain_length,xsec
680 $ , scale_vf(3)
681 common /cvflow_i/ icvflow,iavflow
682 common /cvflow_c/ chv(3)
683 character*1 chv
684c
685 real bd_vflow,dt_vflow
686 save bd_vflow,dt_vflow
687 data bd_vflow,dt_vflow /-99.,-99./
688
689 logical ifcomp
690
691c Check list:
692
693c param (55) -- volume flow rate, if nonzero
694c forcing in X? or in Z?
695
696 ntot1 = lx1*ly1*lz1*nelv
697 ntot2 = lx2*ly2*lz2*nelv
698
699 if (param(55).eq.0.) return
700 if (kx1.eq.1) then
701 write(6,*) 'ABORT. Recompile vol_flow with kx1=lx1, etc.'
702 call exitt
703 endif
704
705 icvflow = 1 ! Default flow dir. = X
706 if (param(54).ne.0) icvflow = abs(param(54))
707 iavflow = 0 ! Determine flow rate
708 if (param(54).lt.0) iavflow = 1 ! from mean velocity
709 flow_rate = param(55)
710
711 chv(1) = 'X'
712 chv(2) = 'Y'
713 chv(3) = 'Z'
714
715c If either dt or the backwards difference coefficient change,
716c then recompute base flow solution corresponding to unit forcing:
717
718 ifcomp = .false.
719 if (dt.ne.dt_vflow.or.bd(1).ne.bd_vflow.or.ifmvbd) ifcomp=.true.
720 if (.not.ifcomp) then
721 ifcomp=.true.
722 do i=1,ntot1
723 if (vdiff (i,1,1,1,1).ne.vdc(i,1)) goto 20
724 if (vtrans(i,1,1,1,1).ne.vdc(i,2)) goto 20
725 enddo
726 ifcomp=.false. ! If here, then vdiff/vtrans unchanged.
727 20 continue
728 endif
729
730 call copy(vdc(1,1),vdiff (1,1,1,1,1),ntot1)
731 call copy(vdc(1,2),vtrans(1,1,1,1,1),ntot1)
732 dt_vflow = dt
733 bd_vflow = bd(1)
734
735 if (ifcomp) call compute_vol_soln_ms(vxc,vyc,vzc,prc)
736
737 if (icvflow.eq.1)
738 $ current_flow=glsc2_ms(vx,bm1ms,ntot1)/domain_length ! for X
739 if (icvflow.eq.2)
740 $ current_flow=glsc2_ms(vy,bm1ms,ntot1)/domain_length ! for Y
741 if (icvflow.eq.3)
742 $ current_flow=glsc2_ms(vz,bm1ms,ntot1)/domain_length ! for Z
743 volvm1ms = glsum_ms(bm1ms,ntot1)
744
745 if (iavflow.eq.1) then
746 xsec = volvm1ms / domain_length
747 flow_rate = param(55)*xsec
748 endif
749
750 delta_flow = flow_rate-current_flow
751
752c Note, this scale factor corresponds to FFX, provided FFX has
753c not also been specified in userf. If ffx is also specified
754c in userf then the true FFX is given by ffx_userf + scale.
755
756 scale = delta_flow/base_flow
757 scale_vf(icvflow) = scale
758 if (nio.eq.0) write(6,1) istep,chv(icvflow)
759 $ ,time,scale,delta_flow,current_flow,flow_rate
760 1 format(i10,' volflow ',a1,11x,1p5e12.4)
761
762 call add2s2(vx,vxc,scale,ntot1)
763 call add2s2(vy,vyc,scale,ntot1)
764 call add2s2(vz,vzc,scale,ntot1)
765 call add2s2(pr,prc,scale,ntot2)
766
767 return
768 end
769c-----------------------------------------------------------------------
770 subroutine compute_vol_soln_ms(vxc,vyc,vzc,prc)
771c
772c Compute the solution to the time-dependent Stokes problem
773c with unit forcing, and find associated flow rate.
774c
775c pff 2/28/98
776c
777 include 'SIZE'
778 include 'TOTAL'
779c
780 real vxc(lx1,ly1,lz1,lelv)
781 $ , vyc(lx1,ly1,lz1,lelv)
782 $ , vzc(lx1,ly1,lz1,lelv)
783 $ , prc(lx2,ly2,lz2,lelv)
784c
785 common /cvflow_r/ flow_rate,base_flow,domain_length,xsec
786 $ , scale_vf(3)
787 common /cvflow_i/ icvflow,iavflow
788 common /cvflow_c/ chv(3)
789 character*1 chv
790c
791 integer icalld
792 save icalld
793 data icalld/0/
794
795c
796c
797 ntot1 = lx1*ly1*lz1*nelv
798 if (icalld.eq.0) then
799 icalld=icalld+1
800 xlmin = glmin_ms(xm1,ntot1)
801 xlmax = glmax_ms(xm1,ntot1)
802 ylmin = glmin_ms(ym1,ntot1) ! for Y!
803 ylmax = glmax_ms(ym1,ntot1)
804 zlmin = glmin_ms(zm1,ntot1) ! for Z!
805 zlmax = glmax_ms(zm1,ntot1)
806c
807 if (icvflow.eq.1) domain_length = xlmax - xlmin
808 if (icvflow.eq.2) domain_length = ylmax - ylmin
809 if (icvflow.eq.3) domain_length = zlmax - zlmin
810 endif
811c
812 if (ifsplit) then
813c call plan2_vol(vxc,vyc,vzc,prc)
814 call plan4_vol_ms(vxc,vyc,vzc,prc)
815 else
816 call plan3_vol_ms(vxc,vyc,vzc,prc)
817 endif
818c
819c Compute base flow rate
820c
821 if (icvflow.eq.1)
822 $ base_flow = glsc2_ms(vxc,bm1ms,ntot1)/domain_length
823 if (icvflow.eq.2)
824 $ base_flow = glsc2_ms(vyc,bm1ms,ntot1)/domain_length
825 if (icvflow.eq.3)
826 $ base_flow = glsc2_ms(vzc,bm1ms,ntot1)/domain_length
827c
828 if (nio.eq.0 .and. loglevel.gt.0) write(6,1)
829 $ istep,chv(icvflow),base_flow,domain_length,flow_rate
830 1 format(i11,' basflow ',a1,11x,1p3e13.4)
831c
832 return
833 end
834c-----------------------------------------------------------------------
835 subroutine plan3_vol_ms(vxc,vyc,vzc,prc)
836c
837c Compute pressure and velocity using fractional step method.
838c (PLAN3).
839c
840c
841 include 'SIZE'
842 include 'TOTAL'
843c
844 real vxc(lx1,ly1,lz1,lelv)
845 $ , vyc(lx1,ly1,lz1,lelv)
846 $ , vzc(lx1,ly1,lz1,lelv)
847 $ , prc(lx2,ly2,lz2,lelv)
848C
849 COMMON /SCRNS/ rw1 (LX1,LY1,LZ1,LELV)
850 $ , rw2 (LX1,LY1,LZ1,LELV)
851 $ , rw3 (LX1,LY1,LZ1,LELV)
852 $ , dv1 (LX1,LY1,LZ1,LELV)
853 $ , dv2 (LX1,LY1,LZ1,LELV)
854 $ , dv3 (LX1,LY1,LZ1,LELV)
855 $ , RESPR (LX2,LY2,LZ2,LELV)
856 COMMON /SCRVH/ H1 (LX1,LY1,LZ1,LELV)
857 $ , H2 (LX1,LY1,LZ1,LELV)
858 COMMON /SCRHI/ H2INV (LX1,LY1,LZ1,LELV)
859 common /cvflow_i/ icvflow,iavflow
860
861 real vxcbc(lx1,ly1,lz1,lelv)
862 real vycbc(lx1,ly1,lz1,lelv)
863 real vzcbc(lx1,ly1,lz1,lelv)
864 REAL vxcp (LX1,LY1,LZ1,LELV)
865 REAL dvxc (LX1,LY1,LZ1,LELV)
866 REAL vycp (LX1,LY1,LZ1,LELV)
867 REAL dvyc (LX1,LY1,LZ1,LELV)
868 REAL vzcp (LX1,LY1,LZ1,LELV)
869 REAL dvzc (LX1,LY1,LZ1,LELV)
870 common /cvflow_nn/ vxcbc,vycbc,vzcbc
871 real resbc(lx1*ly1*lz1*lelv,ldim+1)
872c
873c
874c Compute velocity, 1st part
875c
876 n = lx1*ly1*lz1*nelv
877 ntot1 = lx1*ly1*lz1*nelv
878 ntot2 = lx2*ly2*lz2*nelv
879 ifield = 1
880c
881 call opzero(vxcbc,vycbc,vzcbc)
882 call opzero(vxc,vyc,vzc)
883
884 ngeompv = 10
885 do ictr = 1,ngeompv
886 if (icvflow.eq.1) then
887 call copy (rw1,bm1,ntot1)
888 call rzero (rw2,ntot1)
889 call rzero (rw3,ntot1)
890 elseif (icvflow.eq.2) then
891 call rzero (rw1,ntot1)
892 call copy (rw2,bm1,ntot1)
893 call rzero (rw3,ntot1)
894 else
895 call rzero (rw1,ntot1) ! Z-flow!
896 call rzero (rw2,ntot1) ! Z-flow!
897 call copy (rw3,bm1,ntot1) ! Z-flow!
898 endif
899
900 if (ictr.eq.1) then
901 intype = -1
902 call sethlm (h1,h2,intype)
903 call ophinv (vxc,vyc,vzc,rw1,rw2,rw3,h1,h2,tolhv,nmxh)
904 call ssnormd (vxc,vyc,vzc)
905 else
906 intype = -1
907 call sethlm (h1,h2,intype)
908
909 call copy(vxcp,vxc,lx1*ly1*lz1*nelv)
910 call copy(vycp,vyc,lx1*ly1*lz1*nelv)
911 if (ldim.eq.3) call copy(vzcp,vzc,lx1*ly1*lz1*nelv)
912
913 call neknek_xfer_fld(vxc,vxcbc)
914 call neknek_xfer_fld(vyc,vycbc)
915 if (ldim.eq.3) call neknek_xfer_fld(vzc,vzcbc)
916
917 call ophx(resbc(1,1),resbc(1,2),resbc(1,3),
918 $ vxcbc,vycbc,vzcbc,h1,h2)
919 call sub2(rw1,resbc(1,1),ntot1)
920 call sub2(rw2,resbc(1,2),ntot1)
921 if (ldim.eq.3) call sub2(rw3,resbc(1,3),ntot1)
922 call ophinv (vxc,vyc,vzc,rw1,rw2,rw3,h1,h2,tolhv,nmxh)
923 call opadd2(vxc,vyc,vzc,vxcbc,vycbc,vzcbc)
924 call ssnormd (vxc,vyc,vzc)
925 endif
926c
927c
928c Compute pressure (from "incompr")
929c
930 intype = 1
931 dtinv = 1./dt
932c
933 call rzero (h1,ntot1)
934 call copy (h2,vtrans(1,1,1,1,ifield),ntot1)
935 call cmult (h2,dtinv,ntot1)
936 call invers2 (h2inv,h2,ntot1)
937 call opdiv (respr,vxc,vyc,vzc)
938 call chsign (respr,ntot2)
939 call ortho (respr)
940c
941c
942c Set istep=0 so that h1/h2 will be re-initialized in eprec
943 i_tmp = istep
944 istep = 0
945 call esolver (respr,h1,h2,h2inv,intype)
946 istep = i_tmp
947c
948 call opgradt (rw1,rw2,rw3,respr)
949 call opbinv (dv1,dv2,dv3,rw1,rw2,rw3,h2inv)
950 call opadd2 (vxc,vyc,vzc,dv1,dv2,dv3)
951c
952 call cmult2 (prc,respr,bd(1),ntot2)
953
954 call sub3(dvxc,vxcp,vxc,ntot1)
955 call sub3(dvyc,vycp,vyc,ntot1)
956 call sub3(dvzc,vzcp,vzc,ntot1)
957 dvxmax = glamax_ms(dvxc,ntot1)
958 dvymax = glamax_ms(dvyc,ntot1)
959 dvzmax = glamax_ms(dvzc,ntot1)
960 if (nio.eq.0)
961 $ write(6,'(i2,i8,i4,1p4e13.4,a11)') idsess,istep,ictr,time,
962 $ dvxmax,dvymax,dvzmax,' del-vol-vxy'
963 call neknekgsync()
964 enddo
965
966 call neknekgsync()
967c
968
969 if (istep.eq.3) ifxyo = .true.
970 if (istep.eq.3) call outpost(vxc,vyc,vzc,prc,t,'cor')
971 if (istep.eq.3) ifxyo = .false.
972 return
973 end
974c-----------------------------------------------------------------------
975 subroutine plan4_vol_ms(vxc,vyc,vzc,prc)
976
977c Compute pressure and velocity using fractional step method.
978c (Tombo splitting scheme).
979 include 'SIZE'
980 include 'TOTAL'
981 include 'NEKNEK'
982
983 real vxc(lx1,ly1,lz1,lelv)
984 $ , vyc(lx1,ly1,lz1,lelv)
985 $ , vzc(lx1,ly1,lz1,lelv)
986 $ , prc(lx2,ly2,lz2,lelv)
987
988 common /scrns/ resv1 (lx1,ly1,lz1,lelv)
989 $ , resv2 (lx1,ly1,lz1,lelv)
990 $ , resv3 (lx1,ly1,lz1,lelv)
991 $ , respr (lx2*ly2*lz2,lelv)
992 $ , TA1 (lx1*ly1*lz1*lelv)
993 $ , TA2 (lx1*ly1*lz1*lelv)
994 $ , TA3 (lx1*ly1*lz1*lelv)
995 $ , WA1 (lx1*ly1*lz1*lelv)
996 $ , WA2 (lx1*ly1*lz1*lelv)
997 $ , WA3 (lx1*ly1*lz1*lelv)
998 common /scrvh/ h1 (lx1,ly1,lz1,lelv)
999 $ , h2 (lx1,ly1,lz1,lelv)
1000 COMMON /SCRMG/ W1 (LX1*LY1*LZ1,LELV)
1001 $ , W2 (LX1*LY1*LZ1,LELV)
1002 $ , W3 (LX1*LY1*LZ1,LELV)
1003
1004 common /cvflow_i/ icvflow,iavflow
1005
1006 real vxcbc(lx1,ly1,lz1,lelv)
1007 real vycbc(lx1,ly1,lz1,lelv)
1008 real vzcbc(lx1,ly1,lz1,lelv)
1009 REAL vxcp (LX1,LY1,LZ1,LELV)
1010 REAL dvxc (LX1,LY1,LZ1,LELV)
1011 REAL vycp (LX1,LY1,LZ1,LELV)
1012 REAL dvyc (LX1,LY1,LZ1,LELV)
1013 REAL vzcp (LX1,LY1,LZ1,LELV)
1014 REAL dvzc (LX1,LY1,LZ1,LELV)
1015 common /cvflow_nn/ vxcbc,vycbc,vzcbc
1016 real resbc(lx1*ly1*lz1*lelv,ldim+1)
1017
1018 CHARACTER CB*3
1019
1020
1021 n = lx1*ly1*lz1*nelv
1022 NXYZ1 = lx1*ly1*lz1
1023 NTOT1 = NXYZ1*NELV
1024
1025 ngeompv = 10
1026 do ictr = 1,ngeompv
1027 call invers2 (h1,vtrans,n)
1028 call rzero (h2, n)
1029 if (ictr.eq.1) then
1030 call opzero(vxc,vyc,vzc)
1031 else
1032 call neknek_xfer_fld(vxc,vxcbc)
1033 call neknek_xfer_fld(vyc,vycbc)
1034 if (ldim.eq.3) call neknek_xfer_fld(vzc,vzcbc)
1035 endif
1036
1037c Compute pressure
1038 if (icvflow.eq.1) call cdtp(respr,h1,rxm2,sxm2,txm2,1)
1039 if (icvflow.eq.2) call cdtp(respr,h1,rym2,sym2,tym2,1)
1040 if (icvflow.eq.3) call cdtp(respr,h1,rzm2,szm2,tzm2,1)
1041
1042 dtbd = BD(1)/DT
1043
1044C surface terms
1045 DO 100 IEL=1,NELV
1046 DO 300 IFC=1,2*ldim
1047 CALL RZERO (W1(1,IEL),NXYZ1)
1048 CALL RZERO (W2(1,IEL),NXYZ1)
1049 IF (ldim.EQ.3)
1050 $ CALL RZERO (W3(1,IEL),NXYZ1)
1051 CB = CBC(IFC,IEL,IFIELD)
1052 IF (CB(1:1).EQ.'v'.and.intflag(ifc,iel).eq.1) then
1053 CALL FACCL3
1054 $ (W1(1,IEL),vxcbc(1,1,1,IEL),UNX(1,1,IFC,IEL),IFC)
1055 CALL FACCL3
1056 $ (W2(1,IEL),vycbc(1,1,1,IEL),UNY(1,1,IFC,IEL),IFC)
1057 IF (ldim.EQ.3)
1058 $ CALL FACCL3
1059 $ (W3(1,IEL),vzcbc(1,1,1,IEL),UNZ(1,1,IFC,IEL),IFC)
1060 ENDIF
1061 CALL ADD2 (W1(1,IEL),W2(1,IEL),NXYZ1)
1062 IF (ldim.EQ.3)
1063 $ CALL ADD2 (W1(1,IEL),W3(1,IEL),NXYZ1)
1064 CALL FACCL2 (W1(1,IEL),AREA(1,1,IFC,IEL),IFC)
1065 IF (CB(1:1).EQ.'v'.and.intflag(ifc,iel).eq.1) then
1066 CALL CMULT(W1(1,IEL),dtbd,NXYZ1)
1067 endif
1068 CALL SUB2 (RESPR(1,IEL),W1(1,IEL),NXYZ1)
1069 300 CONTINUE
1070 100 CONTINUE
1071
1072
1073 call ortho (respr)
1074 call ctolspl (tolspl,respr)
1075
1076 call hmholtz ('PRES',prc,respr,h1,h2,pmask,vmult,
1077 $ imesh,tolspl,nmxh,1)
1078 call ortho (prc)
1079
1080C Compute velocity
1081ccccc
1082 call opgrad (resv1,resv2,resv3,prc)
1083 if (ifaxis) call col2 (resv2,omask,n)
1084 call opchsgn (resv1,resv2,resv3)
1085
1086 if (icvflow.eq.1) call add2col2(resv1,v1mask,bm1,n) ! add forcing
1087 if (icvflow.eq.2) call add2col2(resv2,v2mask,bm1,n)
1088 if (icvflow.eq.3) call add2col2(resv3,v3mask,bm1,n)
1089 if (ifexplvis) call split_vis ! split viscosity into exp/imp part
1090
1091 call copy(vxcp,vxc,lx1*ly1*lz1*nelv)
1092 call copy(vycp,vyc,lx1*ly1*lz1*nelv)
1093 call copy(vzcp,vzc,lx1*ly1*lz1*nelv)
1094
1095 intype = -1
1096 call sethlm (h1,h2,intype)
1097 call ophx(resbc(1,1),resbc(1,2),resbc(1,3),
1098 $ vxcbc,vycbc,vzcbc,h1,h2)
1099 call sub2(resv1,resbc(1,1),n)
1100 call sub2(resv2,resbc(1,2),n)
1101 if (ldim.eq.3) call sub2(resv3,resbc(1,3),n)
1102 call ophinv (vxc,vyc,vzc,resv1,resv2,resv3,h1,h2,tolhv,nmxh)
1103 call opadd2(vxc,vyc,vzc,vxcbc,vycbc,vzcbc)
1104
1105 call sub3(dvxc,vxcp,vxc,n)
1106 call sub3(dvyc,vycp,vyc,n)
1107 call sub3(dvzc,vzcp,vzc,n)
1108 dvxmax = glamax_ms(dvxc,n)
1109 dvymax = glamax_ms(dvyc,n)
1110 dvzmax = glamax_ms(dvzc,n)
1111 if (ifexplvis) call redo_split_vis ! restore vdiff
1112 if (nid.eq.0)
1113 $ write(6,'(i2,i8,i4,1p4e13.4,a11)') idsess,istep,ictr,time,
1114 $ dvxmax,dvymax,dvzmax,' del-vol-vxy'
1115
1116 enddo
1117
1118 return
1119 end
1120c-----------------------------------------------------------------------
1121
Note: See TracBrowser for help on using the repository browser.