KTH framework for Nek5000 toolboxes; testing version  0.0.1
multimesh.f
Go to the documentation of this file.
1 c-----------------------------------------------------------------------
2 c Routines for multidomain (neknek) simulations.
3 c
4 c References:
5 c
6 c "A spectrally accurate method for overlapping grid solution of
7 c incompressible Navier–Stokes equations" Brandon E. Merrill,
8 c Yulia T. Peet, Paul F. Fischer, and James W. Lottes, J. Comp. Phys.
9 c 307 (2016) 60-93.
10 c
11 c "Stability analysis of interface temporal discretization in grid
12 c overlapping methods," Y. Peet, P.F. Fischer, SIAM J. Numer. Anal.
13 c 50 (6) (2012) 3375–3401.
14 c-----------------------------------------------------------------------
15  subroutine setup_neknek_wts
16 
17  include 'SIZE'
18  include 'TOTAL'
19 
20  ntot1 = lx1*ly1*lz1*nelt
21 c Initialize unity partition function to 1
22  call rone(upf,ntot1)
23  call col3(bm1ms,bm1,upf,ntot1)
24 
25  return
26  end
27 c-------------------------------------------------------------
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 
68 c Figure out the displacement for the first mesh
69  call setup_int_neknek(dxf,dyf,dzf) !sets up interpolation for 2 meshes
70 
71 c exchange_points finds the processor and element number at
72 c 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
82 c-------------------------------------------------------------
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 
92 c Set interpolation flag: points with boundary condition = 'int'
93 c get intflag=1.
94 c
95 c 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 '
111 c if (cb.eq.'inp') cbc(f,e,ifield)='on ' ! Pressure
112 c 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 
119 c 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
126 c------------------------------------------------------------------------
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 
143 c Order of extrpolation is contolled by the parameter NINTER contained
144 c in NEKNEK. First order interface extrapolation, NINTER=1 (time lagging)
145 c is activated. It is unconditionally stable. If you want to use
146 c higher-order interface extrapolation schemes, you need to increase
147 c 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
173 c---------------------------------------------------------------------
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
213 c-----------------------------------------------------------------------
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
231 c-----------------------------------------------------------------------
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
242 c THIS ROUTINE DISPLACES THE FIRST MESH AND SETUPS THE FINDPTS
243 c THE MESH IS DISPLACED BACK TO ORIGINAL POSITION IN EXCH_POINTS
244 
245 c 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 
251 c 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 
260 c Displace MESH 1
261  ntot = lx1*ly1*lz1*nelt
262  if (idsess.eq.0) then
263  call cadd(xm1,-dxf,ntot)
264  endif
265 
266 c 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
282 c-----------------------------------------------------------------------
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 
302 c Look for boundary points with Diriclet b.c. (candidates for
303 c 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 
316 c Setup arrays of x,y,zs to send to findpts and indices of boundary
317 c 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 
360 c 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 
372 c 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 
381 c 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
416 c-----------------------------------------------------------------------
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 
442 c 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 
453 c Now we can transfer this information to valint array from which
454 c 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
473 c--------------------------------------------------------------------------
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 
482 c 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
492 c--------------------------------------------------------------------------
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
507 C--------------------------------------------------------------------------
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 
518 c Now we can transfer this information to valint array from which
519 c 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
528 C--------------------------------------------------------------------------
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/
540 c assume that this routine is called at the end of bcdirvc
541 c where all the boundary condition data has been read in for
542 c 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
574 c if (nid.eq.0) write(6,104) idsess,istep,time,dqg,aqg,gamma
575 c 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
597 c-----------------------------------------------------------------------
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
629 c-----------------------------------------------------------------------
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
652 c-----------------------------------------------------------------------
653 c-----------------------------------------------------------------------
654  subroutine vol_flow_ms
655 c
656 c
657 c Adust flow volume at end of time step to keep flow rate fixed by
658 c adding an appropriate multiple of the linear solution to the Stokes
659 c problem arising from a unit forcing in the X-direction. This assumes
660 c that the flow rate in the X-direction is to be fixed (as opposed to Y-
661 c or Z-) *and* that the periodic boundary conditions in the X-direction
662 c occur at the extreme left and right ends of the mesh.
663 c
664 c pff 6/28/98
665 c
666  include 'SIZE'
667  include 'TOTAL'
668 c
669 c Swap the comments on these two lines if you don't want to fix the
670 c flow rate for periodic-in-X (or Z) flow problems.
671 c
672  parameter(kx1=lx1,ky1=ly1,kz1=lz1,kx2=lx2,ky2=ly2,kz2=lz2)
673 c
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
684 c
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 
691 c Check list:
692 
693 c param (55) -- volume flow rate, if nonzero
694 c 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 
715 c If either dt or the backwards difference coefficient change,
716 c 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 
752 c Note, this scale factor corresponds to FFX, provided FFX has
753 c not also been specified in userf. If ffx is also specified
754 c 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
769 c-----------------------------------------------------------------------
770  subroutine compute_vol_soln_ms(vxc,vyc,vzc,prc)
771 c
772 c Compute the solution to the time-dependent Stokes problem
773 c with unit forcing, and find associated flow rate.
774 c
775 c pff 2/28/98
776 c
777  include 'SIZE'
778  include 'TOTAL'
779 c
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)
784 c
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
790 c
791  integer icalld
792  save icalld
793  data icalld/0/
794 
795 c
796 c
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)
806 c
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
811 c
812  if (ifsplit) then
813 c 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
818 c
819 c Compute base flow rate
820 c
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
827 c
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)
831 c
832  return
833  end
834 c-----------------------------------------------------------------------
835  subroutine plan3_vol_ms(vxc,vyc,vzc,prc)
836 c
837 c Compute pressure and velocity using fractional step method.
838 c (PLAN3).
839 c
840 c
841  include 'SIZE'
842  include 'TOTAL'
843 c
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)
848 C
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)
872 c
873 c
874 c Compute velocity, 1st part
875 c
876  n = lx1*ly1*lz1*nelv
877  ntot1 = lx1*ly1*lz1*nelv
878  ntot2 = lx2*ly2*lz2*nelv
879  ifield = 1
880 c
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
926 c
927 c
928 c Compute pressure (from "incompr")
929 c
930  intype = 1
931  dtinv = 1./dt
932 c
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)
940 c
941 c
942 c 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
947 c
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)
951 c
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()
967 c
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
974 c-----------------------------------------------------------------------
975  subroutine plan4_vol_ms(vxc,vyc,vzc,prc)
976 
977 c Compute pressure and velocity using fractional step method.
978 c (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 
1037 c 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 
1044 C 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 
1080 C Compute velocity
1081 ccccc
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
1120 c-----------------------------------------------------------------------
1121 
nekasgn
subroutine nekasgn(ix, iy, iz, e)
Definition: bdry.f:1121
exitt
void exitt()
Definition: comm_mpi.f:604
nekgsync
subroutine nekgsync()
Definition: comm_mpi.f:502
exitti
subroutine exitti(stringi, idata)
Definition: comm_mpi.f:535
neknekgsync
subroutine neknekgsync()
Definition: comm_mpi.f:1401
dnekclock
real *8 function dnekclock()
Definition: comm_mpi.f:393
dnekclock_sync
real *8 function dnekclock_sync()
Definition: comm_mpi.f:401
facind
subroutine facind(kx1, kx2, ky1, ky2, kz1, kz2, nx, ny, nz, iface)
Definition: connect1.f:1028
dsset
subroutine dsset(nx, ny, nz)
Definition: connect1.f:553
scale
subroutine scale(xyzl, nl)
Definition: connect2.f:602
hmholtz
subroutine hmholtz(name, u, rhs, h1, h2, mask, mult, imsh, tli, maxit, isd)
Definition: hmholtz.f:3
geom_reset
subroutine geom_reset(icall)
Definition: ic.f:1802
ophinv
subroutine ophinv(o1, o2, o3, i1, i2, i3, h1, h2, tolh, nmxhi)
Definition: induct.f:1023
opzero
subroutine opzero(ux, uy, uz)
Definition: induct.f:406
cadd
subroutine cadd(a, const, n)
Definition: math.f:326
col3
subroutine col3(a, b, c, n)
Definition: math.f:598
sub2
subroutine sub2(a, b, n)
Definition: math.f:164
invers2
subroutine invers2(a, b, n)
Definition: math.f:51
glmin
function glmin(a, n)
Definition: math.f:973
add2col2
subroutine add2col2(a, b, c, n)
Definition: math.f:1565
glmin_ms
function glmin_ms(a, n)
Definition: math.f:1633
iglsum
function iglsum(a, n)
Definition: math.f:926
glsum_ms
function glsum_ms(a, n)
Definition: math.f:1681
add2
subroutine add2(a, b, n)
Definition: math.f:622
col2
subroutine col2(a, b, n)
Definition: math.f:564
add2s2
subroutine add2s2(a, b, c1, n)
Definition: math.f:690
copy
subroutine copy(a, b, n)
Definition: math.f:260
iglmax
function iglmax(a, n)
Definition: math.f:913
izero
subroutine izero(a, n)
Definition: math.f:215
iglmin_ms
function iglmin_ms(a, n)
Definition: math.f:1741
glmax_ms
function glmax_ms(a, n)
Definition: math.f:1657
glsum
function glsum(x, n)
Definition: math.f:861
sub3
subroutine sub3(a, b, c, n)
Definition: math.f:175
cmult
subroutine cmult(a, const, n)
Definition: math.f:315
iglmax_ms
function iglmax_ms(a, n)
Definition: math.f:1729
glamax_ms
function glamax_ms(a, n)
Definition: math.f:1669
rone
subroutine rone(a, n)
Definition: math.f:230
glsc2_ms
function glsc2_ms(a, b, n)
Definition: math.f:1717
rzero
subroutine rzero(a, n)
Definition: math.f:208
glmax
function glmax(a, n)
Definition: math.f:960
chsign
subroutine chsign(a, n)
Definition: math.f:305
compute_vol_soln_ms
subroutine compute_vol_soln_ms(vxc, vyc, vzc, prc)
Definition: multimesh.f:771
plan4_vol_ms
subroutine plan4_vol_ms(vxc, vyc, vzc, prc)
Definition: multimesh.f:976
field_eval
subroutine field_eval(fieldout, fieldstride, fieldin)
Definition: multimesh.f:475
vol_flow_ms
subroutine vol_flow_ms
Definition: multimesh.f:655
neknekmv
subroutine neknekmv()
Definition: multimesh.f:215
exchange_points
subroutine exchange_points(dxf, dyf, dzf)
Definition: multimesh.f:284
plan3_vol_ms
subroutine plan3_vol_ms(vxc, vyc, vzc, prc)
Definition: multimesh.f:836
fix_surface_flux
subroutine fix_surface_flux
Definition: multimesh.f:530
nekneksanchk
subroutine nekneksanchk
Definition: multimesh.f:494
chk_outflow
subroutine chk_outflow
Definition: multimesh.f:175
rescale_x_ms
subroutine rescale_x_ms(x, x0, x1)
Definition: multimesh.f:631
neknek_xfer_fld
subroutine neknek_xfer_fld(u, ui)
Definition: multimesh.f:509
bcopy
subroutine bcopy
Definition: multimesh.f:128
neknek_exchange
subroutine neknek_exchange
Definition: multimesh.f:418
set_intflag
subroutine set_intflag
Definition: multimesh.f:84
surface_flux_area
subroutine surface_flux_area(dq, aq, qx, qy, qz, e, f, w)
Definition: multimesh.f:599
setup_int_neknek
subroutine setup_int_neknek(dxf, dyf, dzf)
Definition: multimesh.f:233
neknek_setup
subroutine neknek_setup
Definition: multimesh.f:29
setup_neknek_wts
subroutine setup_neknek_wts
Definition: multimesh.f:16
esolver
subroutine esolver(RES, H1, H2, H2INV, INTYPE)
Definition: navier0.f:2
opgradt
subroutine opgradt(outx, outy, outz, inpfld)
Definition: navier1.f:4096
cdtp
subroutine cdtp(dtx, x, rm2, sm2, tm2, isd)
Definition: navier1.f:331
ortho
subroutine ortho(respr)
Definition: navier1.f:224
opbinv
subroutine opbinv(out1, out2, out3, inp1, inp2, inp3, h2inv)
Definition: navier1.f:776
opdiv
subroutine opdiv(outfld, inpx, inpy, inpz)
Definition: navier1.f:4065
ophx
subroutine ophx(out1, out2, out3, inp1, inp2, inp3, h1, h2)
Definition: navier1.f:743
ctolspl
subroutine ctolspl(tolspl, respr)
Definition: navier1.f:194
opadd2
subroutine opadd2(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2350
opchsgn
subroutine opchsgn(a, b, c)
Definition: navier1.f:2464
opgrad
subroutine opgrad(out1, out2, out3, inp)
Definition: navier1.f:298
mappr
subroutine mappr(pm1, pm2, pa, pb)
Definition: navier5.f:237
fix_geom
subroutine fix_geom
Definition: navier5.f:2322
split_vis
subroutine split_vis
Definition: plan4.f:459
redo_split_vis
subroutine redo_split_vis
Definition: plan4.f:477
outpost
subroutine outpost(v1, v2, v3, vp, vt, name3)
Definition: prepost.f:1378
ssnormd
subroutine ssnormd(DV1, DV2, DV3)
Definition: ssolv.f:629
faccl2
subroutine faccl2(a, b, iface1)
Definition: subs1.f:909
faddcl3
subroutine faddcl3(a, b, c, iface1)
Definition: subs1.f:978
faccl3
subroutine faccl3(a, b, c, iface1)
Definition: subs1.f:944
sethlm
subroutine sethlm(h1, h2, intloc)
Definition: subs1.f:1014
cmult2
subroutine cmult2(A, B, CONST, N)
Definition: subs2.f:341