[Nek5000-users] [Re] Userdat2 and BC

nek5000-users at lists.mcs.anl.gov nek5000-users at lists.mcs.anl.gov
Thu Jun 7 11:15:08 CDT 2018


Thanks Paul, helps a lot!

Juan Pablo

El jue., 7 de jun. de 2018 07:59, <nek5000-users-request at lists.mcs.anl.gov>
escribió:

> Send Nek5000-users mailing list submissions to
>         nek5000-users at lists.mcs.anl.gov
>
> To subscribe or unsubscribe via the World Wide Web, visit
>         https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
> or, via email, send a message with subject or body 'help' to
>         nek5000-users-request at lists.mcs.anl.gov
>
> You can reach the person managing the list at
>         nek5000-users-owner at lists.mcs.anl.gov
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of Nek5000-users digest..."
>
>
> Today's Topics:
>
>    1. Userdat2 and BC (nek5000-users at lists.mcs.anl.gov)
>    2. Re: Userdat2 and BC (nek5000-users at lists.mcs.anl.gov)
>    3. Problem size requires AMG solver (nek5000-users at lists.mcs.anl.gov)
>    4. Initial Condition & Recycling (nek5000-users at lists.mcs.anl.gov)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Wed, 6 Jun 2018 18:13:53 -0400
> From: nek5000-users at lists.mcs.anl.gov
> To: nek5000-users at lists.mcs.anl.gov
> Subject: [Nek5000-users] Userdat2 and BC
> Message-ID:
>         <mailman.3469.1528323237.86639.nek5000-users at lists.mcs.anl.gov>
> Content-Type: text/plain; charset="utf-8"
>
> Dear Neks,
>
> I have the following question: the boundary conditions specified in the usr
> file are applied before or after the mesh deformation in userdat2?
>
> Thanks in advance,
>
> Juan Pablo.
> -------------- next part --------------
> An HTML attachment was scrubbed...
> URL: <
> http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20180606/225ad344/attachment-0001.html
> >
>
> ------------------------------
>
> Message: 2
> Date: Thu, 7 Jun 2018 03:28:16 +0000
> From: nek5000-users at lists.mcs.anl.gov
> To: "nek5000-users at lists.mcs.anl.gov"
>         <nek5000-users at lists.mcs.anl.gov>
> Subject: Re: [Nek5000-users] Userdat2 and BC
> Message-ID:
>         <mailman.3477.1528342098.86639.nek5000-users at lists.mcs.anl.gov>
> Content-Type: text/plain; charset="iso-8859-1"
>
>
> Dear Juan Pablo,
>
>
> Boundary conditions are not set until the calculation starts
>
> (e.g., T, vx, vy, vz values are not set before usrdat2 is called).
>
>
> The BC types are known by the time usrdat2 is called (i.e.,
>
> cbc(f,e,ifield) is defined).
>
>
> Does this help?
>
>
> Paul
>
>
> ________________________________
> From: Nek5000-users <nek5000-users-bounces at lists.mcs.anl.gov> on behalf
> of nek5000-users at lists.mcs.anl.gov <nek5000-users at lists.mcs.anl.gov>
> Sent: Wednesday, June 6, 2018 5:13:53 PM
> To: nek5000-users at lists.mcs.anl.gov
> Subject: [Nek5000-users] Userdat2 and BC
>
> Dear Neks,
>
> I have the following question: the boundary conditions specified in the
> usr file are applied before or after the mesh deformation in userdat2?
>
> Thanks in advance,
>
> Juan Pablo.
>
> -------------- next part --------------
> An HTML attachment was scrubbed...
> URL: <
> http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20180607/5df0344f/attachment-0001.html
> >
>
> ------------------------------
>
> Message: 3
> Date: Thu, 7 Jun 2018 07:39:20 +0000
> From: nek5000-users at lists.mcs.anl.gov
> To: "nek5000-users at lists.mcs.anl.gov"
>         <nek5000-users at lists.mcs.anl.gov>
> Subject: [Nek5000-users] Problem size requires AMG solver
> Message-ID:
>         <mailman.3502.1528372770.86639.nek5000-users at lists.mcs.anl.gov>
> Content-Type: text/plain; charset="us-ascii"
>
> Hi everyone,
>
> I need to run a case in Nek with a large number of elements (~ 1 million).
>
> The same problem with a small number of elements works fine, but on the
> larger grid
> the solver exits with the following error:
>
>     EXIT: Problem size requires AMG solver          1
>
> Going through the code
>
> # if (imode.eq.0 .and. nelgt.gt.350000) call exitti(
> #   $ 'Problem size requires AMG solver$',1)
>
> It seems that over 350000 elements, you need to use the AMG solver.
> Is that correct? Or is there any workaround to not use the AMG.
>
> Kind regards,
> Dante De Santis
>
> Research Consultant Computational Physics 4 Solutions
> [NRG - Research and Innovation Unit]
>
> [Logo_NRG_PMS_355-186_C]
> Westerduinweg 3, 1755 LE  PETTEN
> P.O. Box 25, 1755 ZG  PETTEN
> THE NETHERLANDS
> phone: +31(0) 224 56 4656
> e-mail: desantis at nrg.eu<mailto:desantis at nrg.eu>
> Visit NRG at www.nrg.eu<http://www.nrg.eu/>
>
> -------------- next part --------------
> An HTML attachment was scrubbed...
> URL: <
> http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20180607/4e5b7e29/attachment.html
> >
> -------------- next part --------------
> A non-text attachment was scrubbed...
> Name: image001.jpg
> Type: image/jpeg
> Size: 4560 bytes
> Desc: image001.jpg
> URL: <
> http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20180607/4e5b7e29/attachment.jpg
> >
>
> ------------------------------
>
> Message: 4
> Date: Thu, 7 Jun 2018 10:51:09 +0200
> From: nek5000-users at lists.mcs.anl.gov
> To: nek5000-users at lists.mcs.anl.gov
> Subject: [Nek5000-users] Initial Condition & Recycling
> Message-ID:
>         <mailman.3503.1528372771.86639.nek5000-users at lists.mcs.anl.gov>
> Content-Type: text/plain; charset="utf-8"; Format="flowed"
>
> Dear All,
>
>
> In an attempt to simulate an established turbulent flow over a flat
> plate with a recycling method for the inlet, it has to my attention that
> a (1/7)th Power Law would be suitable as an initial condition. However
> the solver divergences and I think I know why.
>
> Here the layout of the recycling procedure :
>
> 1. call set_inflow_fpt in userchk
>
> 2. If it is the 1st call then call setp_inflow_fpt_setup
>
> 3. call rescale_fpt (which rescales the copied field to a user specified
> mean velocity ubar)
>
> 4.In the rescale_fpt routine we call get_flux_and_area which gives the
> flux and area of a face with cbc == 'v   ' so a user specified velocity
> boundary condition.
>
> 4. a) As the initial condition is wrong once that routine sums the flux
> on the boundary it returns zero in this case and the variable scale
> which is divided by the latter returns infinity !
>
>
> However I set the initial conditions as in many other examples, such as
> the turbInflow. Would have an insight that could lead me to a solution ?
>
>
> Please find attached the related files.
>
>
> Sincerely,
>
> Armand, ONERA - The French Aerospace Lab.
>
>
> -------------- next part --------------
> A non-text attachment was scrubbed...
> Name: amg.dat
> Type: application/octet-stream
> Size: 80000 bytes
> Desc: not available
> URL: <
> http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20180607/89aba907/attachment.obj
> >
> -------------- next part --------------
> A non-text attachment was scrubbed...
> Name: amg_Aff.dat
> Type: application/octet-stream
> Size: 943704 bytes
> Desc: not available
> URL: <
> http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20180607/89aba907/attachment-0001.obj
> >
> -------------- next part --------------
> A non-text attachment was scrubbed...
> Name: amg_AfP.dat
> Type: application/octet-stream
> Size: 1794968 bytes
> Desc: not available
> URL: <
> http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20180607/89aba907/attachment-0002.obj
> >
> -------------- next part --------------
> A non-text attachment was scrubbed...
> Name: amg_W.dat
> Type: application/octet-stream
> Size: 673416 bytes
> Desc: not available
> URL: <
> http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20180607/89aba907/attachment-0003.obj
> >
> -------------- next part --------------
> c
> c     Include file to dimension static arrays
> c     and to set some hardwired run-time parameters
> c
>       integer ldim,lx1,lxd,lx2,lx1m,lelg,lelt,lpmin,lpmax,ldimt
>       integer lpelt,lbelt,toteq,lcvelt
>       integer lelx,lely,lelz,mxprev,lgmres,lorder,lhis
>       integer maxobj,lpert,nsessmax,lxo
>       integer lfdm, ldimt_proj
>
>       ! BASIC
>       parameter (ldim=3)               ! domain dimension (2 or 3)
>       parameter (lx1=8)                ! p-order (avoid uneven and values
> <6)
>       parameter (lxd=8)                ! p-order for over-integration
> (dealiasing)
>       parameter (lx2=lx1-0)            ! p-order for pressure (lx1 or
> lx1-2)
>
>       parameter (lelg=1540)            ! max total number of elements
>       parameter (lpmin=4)              ! min MPI ranks
>       parameter (lpmax=64)             ! max MPI ranks
>       parameter (ldimt=2)              ! max auxiliary fields (temperature
> + scalars)
>
>       ! OPTIONAL
>       parameter (ldimt_proj=1)         ! max auxiliary fields residual
> projection
>       parameter (lhis=100)             ! max history points
>       parameter (maxobj=4)             ! max number of objects
>       parameter (lpert=1)              ! max perturbation modes
>       parameter (toteq=1)              ! max number of conserved scalars
> in cmt
>       parameter (nsessmax=1)           ! max sessions
>       parameter (lxo=lx1)              ! max grid size on output (lxo>=lx1)
>       parameter (lelx=16,lely=12,lelz=8) ! global tensor mesh dimensions
>       parameter (mxprev=30,lgmres=20)  ! max dim of projection & Krylov
> space
>       parameter (lorder=2)             ! max order in time
>
>       parameter (lelt=lelg/lpmin + 2)  ! max number of local elements per
> MPI rank
>       parameter (lx1m=1)               ! polynomial order mesh solver
>       parameter (lbelt=1)              ! mhd
>       parameter (lpelt=1)              ! linear stability
>       parameter (lcvelt=1)             ! cvode
>       parameter (lfdm=0)               ! == 1 for fast diagonalization
> method
>
>       ! INTERNALS
>       include 'SIZE.inc'
> -------------- next part --------------
> A non-text attachment was scrubbed...
> Name: turbLOOT.box
> Type: application/vnd.previewsystems.box
> Size: 792 bytes
> Desc: not available
> URL: <
> http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20180607/89aba907/attachment.box
> >
> -------------- next part --------------
> #
> # nek parameter file
> #
> [GENERAL]
> #startFrom = restart.fld
> stopAt = endTime
> endTime = 50
> #numSteps = 0
>
> dt = 0
> timeStepper = bdf2
> extrapolation = OIFS
> variableDt = yes
> targetCFL = 3.5
>
> writeControl = timeStep #Was runTime
> writeInterval = 1
>
> userParam01 = 25.0    # start time collecting statistics default 100
> userParam02 = 5.0     # writeInterval 1D statistics 20
>
> #dealiasing = no
>
> filtering = hpfrt # set to none in case of Smagorinski
> filterWeight = 10
> filterCutoffRatio = 0.9
>
> [PROBLEMTYPE]
> variableProperties = none # set to yes in case of Smagorinski
> equation = incompNS
>
> [PRESSURE]
> preconditioner = semg_amg
> residualTol = 1e-04
> residualProj = yes
>
> [VELOCITY]
> residualTol = 1e-07
> density = 1
> viscosity = -840000
> -------------- next part --------------
> c-----------------------------------------------------------------------
> c   CASE PARAMETERS:
>
> c start time for averaging
> #define tSTATSTART uparam(1)
>
> c output frequency for statistics
> #define tSTATFREQ  uparam(2)
>
> c mesh dimensions
> #define PI (4.*atan(1.))
> #define NUMBER_ELEMENTS_X 16
> #define NUMBER_ELEMENTS_Y 12
> #define NUMBER_ELEMENTS_Z 8
>
> c-----------------------------------------------------------------------
>       subroutine uservp (ix,iy,iz,ieg)
>       include 'SIZE'
>       include 'TOTAL'
>       include 'NEKUSE'
>
>       common /cdsmag/ ediff(lx1,ly1,lz1,lelv)
>
>       ie     = gllel(ieg)
>       udiff  = ediff(ix,iy,iz,ie)
>       utrans = 1.
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine userf  (ix,iy,iz,ieg)
>       include 'SIZE'
>       include 'TOTAL'
>       include 'NEKUSE'
>
>       ffx = 0.0 ! This value determined from fixed U_b=1 case.
>       ffy = 0.0
>       ffz = 0.0
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine userq  (ix,iy,iz,ieg)
>       include 'SIZE'
>       include 'TOTAL'
>       include 'NEKUSE'
>
>       qvol   = 0.0
>       source = 0.0
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine userchk
>       include 'SIZE'
>       include 'TOTAL'
>       include 'ZPER'  ! for nelx,nely,nelz
>
>       real x0(3)
>       save x0
>
>       integer icalld
>       save    icalld
>       data    icalld /0/
>
>       real atime,timel
>       save atime,timel
>
>       integer ntdump
>       save    ntdump
>
>       common /cdsmag/ ediff(lx1,ly1,lz1,lelv)
>
>       common /plane/  uavg_pl(ly1*lely)
>      $             ,  urms_pl(ly1*lely)
>      $             ,  vrms_pl(ly1*lely)
>      $             ,  wrms_pl(ly1*lely)
>      $             ,  uvms_pl(ly1*lely)
>      $             ,  yy(ly1*lely)
>      $             ,  w1(ly1*lely),w2(ly1*lely)
>
>       common /avg/    uavg(lx1,ly1,lz1,lelv)
>      &             ,  urms(lx1,ly1,lz1,lelv)
>      &             ,  vrms(lx1,ly1,lz1,lelv)
>      &             ,  wrms(lx1,ly1,lz1,lelv)
>      &             ,  uvms(lx1,ly1,lz1,lelv)
>
>       common /ctorq/ dragx(0:maxobj),dragpx(0:maxobj),dragvx(0:maxobj)
>      $             , dragy(0:maxobj),dragpy(0:maxobj),dragvy(0:maxobj)
>      $             , dragz(0:maxobj),dragpz(0:maxobj),dragvz(0:maxobj)
> c
>      $             , torqx(0:maxobj),torqpx(0:maxobj),torqvx(0:maxobj)
>      $             , torqy(0:maxobj),torqpy(0:maxobj),torqvy(0:maxobj)
>      $             , torqz(0:maxobj),torqpz(0:maxobj),torqvz(0:maxobj)
> c
>      $             , dpdx_mean,dpdy_mean,dpdz_mean
>      $             , dgtq(3,4)
>
>       integer e
>       logical ifverbose
>       common /gaaa/    wo1(lx1,ly1,lz1,lelv)
>      &              ,  wo2(lx1,ly1,lz1,lelv)
>      &              ,  wo3(lx1,ly1,lz1,lelv)
>
>       parameter (lt=lx1*ly1*lz1*lelt)
>       common /myoutflow / d(lt),w1n(lt)
>       real dx,dy,dz,ubar
>
>       real tplus
>
>       n=nx1*ny1*nz1*nelv
>       ! do some checks
>       if (istep.eq.0) then
>          if(mod(nely,2).ne.0) then
>            if(nid.eq.0) write(6,*) 'ABORT: nely has to be even!'
>            call exitt
>          endif
>          if(nelx.gt.lelx .or. nely.gt.lely .or. nelz.gt.lelz) then
>            if(nid.eq.0) write(6,*) 'ABORT: nel_xyz > lel_xyz!'
>            call exitt
>          endif
>
>          call set_obj                   ! objects for surface integrals
>          call rzero(x0,3)               ! torque w.r.t. x0
>       endif
>
>       ! Compute eddy viscosity using dynamic smagorinsky model
>       if(ifuservp) then
>         if(nid.eq.0) write(6,*) 'Calculating eddy visosity'
>         do e=1,nelv
>            call eddy_visc(ediff,e)
>         enddo
>         call copy(t,ediff,n)
>       endif
>
>
>       ub = glsc2(vx,bm1,n)/volvm1
>       e2   = glsc3(vy,bm1,vy,n)+glsc3(vz,bm1,vz,n)
>       e2   = e2/volvm1
>       if(nid.eq.0) write(6,2) time,ub,e2
>    2               format(1p3e13.4,' monitor')
>
>       dx=5
>       dy=0.
>       dz=0.
>       ubar = 1.0
>       call set_inflow_fpt(dx,dy,dz,ubar)
>
> c
> c     Below is just for postprocessing ...
> c
>       if (time.lt.tSTATSTART) goto 999 ! don't start averaging
>
>       nelx  = NUMBER_ELEMENTS_X
>       nely  = NUMBER_ELEMENTS_Y
>       nelz  = NUMBER_ELEMENTS_Z
>
>       rho    = 1.
>       dnu    = param(2)
>       A_w    = XLEN * ZLEN
>
>       if(ifoutfld) then
>         if(ldimt.ge.2) call lambda2(t(1,1,1,1,2))
>         if(ldimt.ge.5) call comp_vort3(t(1,1,1,1,3),wo1,wo2,vx,vy,vz)
>       endif
>
>       call torque_calc(1.0,x0,.false.,.false.) ! wall shear
>
>       if(icalld.eq.0) then
>         call rzero(uavg,n)
>         call rzero(urms,n)
>         call rzero(vrms,n)
>         call rzero(wrms,n)
>         call rzero(uvms,n)
>         dragx_avg = 0
>         atime = 0
>         timel = time
>         call planar_average_s(yy,ym1,w1,w2)
>         icalld = 1
>         ntdump = int(time/tSTATFREQ)
>         if(nid.eq.0) write(6,*) 'Start collecting statistics ...'
>       endif
>
>       dtime = time - timel
>       atime = atime + dtime
>
>       if (atime.ne.0. .and. dtime.ne.0.) then
>         beta      = dtime/atime
>         alpha     = 1.-beta
>         ifverbose = .false.
>
>         call avg1(uavg,vx   ,alpha,beta,n,'uavg',ifverbose)
>         call avg2(urms,vx   ,alpha,beta,n,'urms',ifverbose)
>         call avg2(vrms,vy   ,alpha,beta,n,'vrms',ifverbose)
>         call avg2(wrms,vz   ,alpha,beta,n,'wrms',ifverbose)
>         call avg3(uvms,vx,vy,alpha,beta,n,'uvmm',ifverbose)
>
>         dragx_avg = alpha*dragx_avg + beta*0.5*(dragx(1)+dragx(2))
>
>         ! averaging over statistical homogeneous directions (r-t)
>         call planar_average_s(uavg_pl,uavg,w1,w2)
>         call planar_average_s(urms_pl,urms,w1,w2)
>         call planar_average_s(vrms_pl,vrms,w1,w2)
>         call planar_average_s(wrms_pl,wrms,w1,w2)
>         call planar_average_s(uvms_pl,uvms,w1,w2)
>
>         ! average over the domain height
>         m = ny1*nely
>         do i=1,ny1*nely
>            uavg_pl(i) = 0.5 * (uavg_pl(i) + uavg_pl(m-i+1))
>            urms_pl(i) = 0.5 * (urms_pl(i) + urms_pl(m-i+1))
>            vrms_pl(i) = 0.5 * (vrms_pl(i) + vrms_pl(m-i+1))
>            wrms_pl(i) = 0.5 * (wrms_pl(i) + wrms_pl(m-i+1))
>         enddo
>       endif
>
>       ! write statistics to file
>       if (nid.eq.0 .and. istep.gt.0 .and.
>      &   time.gt.(ntdump+1)*tSTATFREQ) then
>
>          tw     = dragx_avg/A_w
>          u_tau  = sqrt(tw/rho)
>          Re_tau = u_tau*DELTA/dnu
>          tplus  = time * Re_tau * u_tau
>
>          ntdump = ntdump + 1
>          write(6,*) 'Dumping statistics ...', tplus, Re_tau
>
>          open(unit=56,file='vel_fluc_prof.dat')
>          write(56,'(A,1pe14.7)') '#time = ', time
>          write(56,'(A)')
>      &   '#  y     y+     uu     vv     ww     uv'
>          open(unit=57,file='mean_prof.dat')
>          write(57,'(A,1pe14.7)') '#time = ', time
>          write(57,'(A)')
>      &   '#  y     y+    Umean'
>
>          m = ny1*nely
>          do i=1,m
>             write(56,3) yy(i)+1
>      &                ,(yy(i)+1)*Re_tau
>      &                , (urms_pl(i)-(uavg_pl(i))**2)/u_tau**2
>      &                , vrms_pl(i)/u_tau**2
>      &                , wrms_pl(i)/u_tau**2
>      &                , uvms_pl(i)/u_tau**2
>             write(57,3)  yy(i) + 1.
>      &                , (yy(i)+1.)*Re_tau
>      &                , uavg_pl(i)/u_tau
>     3       format(1p15e17.9)
>         enddo
>         close(56)
>         close(57)
>       endif
>
>       timel = time
>  999  continue
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine userbc (ix,iy,iz,iside,ieg)
>       include 'SIZE'
>       include 'TOTAL'
>       include 'NEKUSE'
>
>       common /cvelbc/ uin(lx1,ly1,lz1,lelv)
>      $              , vin(lx1,ly1,lz1,lelv)
>      $              , win(lx1,ly1,lz1,lelv)
>
>       integer e,f,eg
>       e = gllel(eg)
>
>       ux=uin(i,j,k,e)
>       uy=vin(i,j,k,e)
>       uz=win(i,j,k,e)
>
>       temp=0.0
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine useric (ix,iy,iz,ieg)
>       include 'SIZE'
>       include 'TOTAL'
>       include 'NEKUSE'
>
>
>       integer idum
>       save    idum
>       data    idum  / 999 /
>
> c     1/7th power law velocity profile for an established boundary layer
>       real Uinf = 8.4
>       real nu = 1e-5
>       real xin = 95.
>       real ymax = 4.
>
>       delta = 0.37*(nu/Uinf)**(1/5.)*(x+xin)**(4/5.)
>
>       ybar = y/ymax
>       ux = Uinf*(ybar/delta)**(1/7.)
>
>       temp=0
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine usrdat   ! This routine to modify element vertices
>       include 'SIZE'      ! _before_ mesh is generated, which
>       include 'TOTAL'     ! guarantees GLL mapping of mesh.
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine usrdat2   ! This routine to modify mesh coordinates
>       include 'SIZE'
>       include 'TOTAL'
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine usrdat3
>       include 'SIZE'
>       include 'TOTAL'
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine set_obj  ! define objects for surface integrals
> c
>       include 'SIZE'
>       include 'TOTAL'
> c
>       integer e,f
> c
> c     Define new objects
> c
>       nobj = 2                  ! for Periodic
>       iobj = 0
>       do ii=nhis+1,nhis+nobj
>          iobj = iobj+1
>          hcode(10,ii) = 'I'
>          hcode( 1,ii) = 'F' ! 'F'
>          hcode( 2,ii) = 'F' ! 'F'
>          hcode( 3,ii) = 'F' ! 'F'
>          lochis(1,ii) = iobj
>       enddo
>       nhis = nhis + nobj
> c
>       if (maxobj.lt.nobj) write(6,*) 'increase maxobj in SIZEu. rm *.o'
>       if (maxobj.lt.nobj) call exitt
> c
>       nxyz = nx1*ny1*nz1
>       do e=1,nelv
>       do f=1,2*ndim
>          if (cbc(f,e,1).eq.'W  ') then
>             iobj = 0
>             if (f.eq.1) iobj=1  ! lower wall
>             if (f.eq.3) iobj=2  ! upper wall
>             if (iobj.gt.0) then
>                nmember(iobj) = nmember(iobj) + 1
>                mem = nmember(iobj)
>                ieg = lglel(e)
>                object(iobj,mem,1) = ieg
>                object(iobj,mem,2) = f
> c              write(6,1) iobj,mem,f,ieg,e,nid,' OBJ'
>     1          format(6i9,a4)
>             endif
> c
>          endif
>       enddo
>       enddo
> c     write(6,*) 'number',(nmember(k),k=1,4)
> c
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine comp_lij(lij,u,v,w,fu,fv,fw,fh,fht,e)
> c
> c     Compute Lij for dynamic Smagorinsky model:
> c                    _   _      _______
> c          L_ij  :=  u_i u_j  - u_i u_j
> c
>       include 'SIZE'
> c
>       integer e
> c
>       real lij(lx1*ly1*lz1,3*ldim-3)
>       real u  (lx1*ly1*lz1,lelv)
>       real v  (lx1*ly1*lz1,lelv)
>       real w  (lx1*ly1*lz1,lelv)
>       real fu (1) , fv (1) , fw (1)
>      $   , fh (1) , fht(1)
>
>       call tens3d1(fu,u(1,e),fh,fht,nx1,nx1)  ! fh x fh x fh x u
>       call tens3d1(fv,v(1,e),fh,fht,nx1,nx1)
>       call tens3d1(fw,w(1,e),fh,fht,nx1,nx1)
>
>       n = nx1*ny1*nz1
>       do i=1,n
>          lij(i,1) = fu(i)*fu(i)
>          lij(i,2) = fv(i)*fv(i)
>          lij(i,3) = fw(i)*fw(i)
>          lij(i,4) = fu(i)*fv(i)
>          lij(i,5) = fv(i)*fw(i)
>          lij(i,6) = fw(i)*fu(i)
>       enddo
>
>       call col3   (fu,u(1,e),u(1,e),n)    !  _______
>       call tens3d1(fv,fu,fh,fht,nx1,nx1)  !  u_1 u_1
>       call sub2   (lij(1,1),fv,n)
>
>       call col3   (fu,v(1,e),v(1,e),n)    !  _______
>       call tens3d1(fv,fu,fh,fht,nx1,nx1)  !  u_2 u_2
>       call sub2   (lij(1,2),fv,n)
>
>       call col3   (fu,w(1,e),w(1,e),n)    !  _______
>       call tens3d1(fv,fu,fh,fht,nx1,nx1)  !  u_3 u_3
>       call sub2   (lij(1,3),fv,n)
>
>       call col3   (fu,u(1,e),v(1,e),n)    !  _______
>       call tens3d1(fv,fu,fh,fht,nx1,nx1)  !  u_1 u_2
>       call sub2   (lij(1,4),fv,n)
>
>       call col3   (fu,v(1,e),w(1,e),n)    !  _______
>       call tens3d1(fv,fu,fh,fht,nx1,nx1)   !  u_2 u_3
>       call sub2   (lij(1,5),fv,n)
>
>       call col3   (fu,w(1,e),u(1,e),n)    !  _______
>       call tens3d1(fv,fu,fh,fht,nx1,nx1)  !  u_3 u_1
>       call sub2   (lij(1,6),fv,n)
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine comp_mij(mij,sij,dg2,fs,fi,fh,fht,nt,e)
> c
> c     Compute Mij for dynamic Smagorinsky model:
> c
> c                     2 _  ____     _______
> c          M_ij  :=  a  S  S_ij  -  S  S_ij
> c
>       include 'SIZE'
> c
>       integer e
> c
>       real mij(lx1*ly1*lz1,3*ldim-3)
>       real dg2(lx1*ly1*lz1,lelv)
>       real fs (1) , fi (1) , fh (1) , fht(1)
>
>       real magS(lx1*ly1*lz1)
>       real sij (lx1*ly1*lz1*ldim*ldim)
>
>       integer imap(6)
>       data imap / 0,4,8,1,5,2 /
>
>       n = nx1*ny1*nz1
>
>       call mag_tensor_e(magS,sij)
>       call cmult(magS,2.0,n)
>
> c     Filter S
>       call tens3d1(fs,magS,fh,fht,nx1,nx1)  ! fh x fh x fh x |S|
>
> c     a2 is the test- to grid-filter ratio, squared
>
>       a2 = nx1-1       ! nx1-1 is number of spaces in grid
>       a2 = a2 /(nt-1)  ! nt-1 is number of spaces in filtered grid
>
>       do k=1,6
>          jj = n*imap(k) + 1
>          call col3   (fi,magS,sij(jj),n)
>          call tens3d1(mij(1,k),fi,fh,fht,nx1,nx1)  ! fh x fh x fh x (|S|
> S_ij)
>          call tens3d1(fi,sij(jj),fh,fht,nx1,nx1)  ! fh x fh x fh x S_ij
>          do i=1,n
>             mij(i,k) = (a2**2 * fs(i)*fi(i) - mij(i,k))*dg2(i,e)
>          enddo
>       enddo
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine eddy_visc(ediff,e)
> c
> c     Compute eddy viscosity using dynamic smagorinsky model
> c
>       include 'SIZE'
>       include 'TOTAL'
>       include 'ZPER'
>
>       real ediff(nx1*ny1*nz1,nelv)
>       integer e
>
>       common /dynsmg/ sij (lx1*ly1*lz1,ldim,ldim)
>      $              , mij (lx1*ly1*lz1,3*ldim-3)
>      $              , lij (lx1*ly1*lz1,3*ldim-3)
>      $              , dg2 (lx1*ly1*lz1,lelv)
>      $              , num (lx1*ly1*lz1,lelv)
>      $              , den (lx1*ly1*lz1,lelv)
>      $              , snrm(lx1*ly1*lz1,lelv)
>      $              , numy(ly1*lely),deny(ly1*lely),yy(ly1*lely)
>       real sij,mij,lij,dg2,num,den,snrm,numy,deny,yy
>
>       parameter(lxyz=lx1*ly1*lz1)
>       common /xzmp0/ ur (lxyz) , us (lxyz) , ut (lxyz)
>       real           vr (lxyz) , vs (lxyz) , vt (lxyz)
>      $     ,         wr (lxyz) , ws (lxyz) , wt (lxyz)
>       common /xzmp1/ w1(lx1*lelv),w2(lx1*lelv)
>
>       !! NOTE CAREFUL USE OF EQUIVALENCE HERE !!
>       equivalence (vr,lij(1,1)),(vs,lij(1,2)),(vt,lij(1,3))
>      $          , (wr,lij(1,4)),(ws,lij(1,5)),(wt,lij(1,6))
>
>       common /sgsflt/ fh(lx1*lx1),fht(lx1*lx1),diag(lx1)
>
>       integer nt
>       save    nt
>       data    nt / -9 /
>
>       ntot = nx1*ny1*nz1
>
>       if (nt.lt.0) call
>      $   set_ds_filt(fh,fht,nt,diag,nx1)! dyn. Smagorinsky filter
>
>       call comp_gije(sij,vx(1,1,1,e),vy(1,1,1,e),vz(1,1,1,e),e)
>       call comp_sije(sij)
>
>       call mag_tensor_e(snrm(1,e),sij)
>       call cmult(snrm(1,e),2.0,ntot)
>
>       call set_grid_spacing(dg2)
>       call comp_mij   (mij,sij,dg2,ur,us,fh,fht,nt,e)
>
>       call comp_lij   (lij,vx,vy,vz,ur,us,ut,fh,fht,e)
>
> c     Compute numerator (ur) & denominator (us) for Lilly contraction
>
>       n = nx1*ny1*nz1
>       do i=1,n
>          ur(i) = mij(i,1)*lij(i,1)+mij(i,2)*lij(i,2)+mij(i,3)*lij(i,3)
>      $      + 2*(mij(i,4)*lij(i,4)+mij(i,5)*lij(i,5)+mij(i,6)*lij(i,6))
>          us(i) = mij(i,1)*mij(i,1)+mij(i,2)*mij(i,2)+mij(i,3)*mij(i,3)
>      $      + 2*(mij(i,4)*mij(i,4)+mij(i,5)*mij(i,5)+mij(i,6)*mij(i,6))
>       enddo
>
> c     smoothing numerator and denominator in time
>       call copy (vr,ur,nx1*nx1*nx1)
>       call copy (vs,us,nx1*nx1*nx1)
>
>       beta1 = 0.0                   ! Temporal averaging coefficients
>       if (istep.gt.1) beta1 = 0.9   ! Retain 90 percent of past
>       beta2 = 1. - beta1
>
>       do i=1,n
>          num (i,e) = beta1*num(i,e) + beta2*vr(i)
>          den (i,e) = beta1*den(i,e) + beta2*vs(i)
>       enddo
>
>
>       if (e.eq.nelv) then  ! planar avg and define nu_tau
>
>          call dsavg(num)   ! average across element boundaries
>          call dsavg(den)
>
>          call planar_average_s      (numy,num,w1,w2)
> c        call wall_normal_average_s (numy,ny1,nely,w1,w2)
>          call planar_fill_s         (num,numy)
>
>          call planar_average_s      (deny,den,w1,w2)
> c        call wall_normal_average_s (deny,ny1,nely,w1,w2)
>          call planar_fill_s         (den,deny)
>
>          call planar_average_s(yy,ym1,w1,w2)
>
> c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> c DIAGNOSTICS ONLY
> c         if (nid.eq.0.and.istep.eq.0) open(unit=55,file='z.z')
> c         if (nid.eq.0.and.mod(istep,10).eq.0) write(55,1)
> c    1    format(/)
> c
> c         ny = ny1*nely
> c         do i=1,ny
> c            cdyn = 0
> c            if (deny(i).gt.0) cdyn = 0.5*numy(i)/deny(i)
> c            cdyn0 = max(cdyn,0.)
> c            if (nid.eq.0.and.mod(istep,10).eq.0) write(55,6)
> c     $         istep,i,time,yy(i),cdyn0,cdyn,numy(i),deny(i)
> c    6       format(i6,i4,1p6e12.4)
> c         enddo
> c DIAGNOSTICS ONLY
> c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
>
>          ntot = nx1*ny1*nz1*nelv
>          do i=1,ntot
>             cdyn = 0
>             if (den(i,1).gt.0) cdyn = 0.5*num(i,1)/den(i,1)
>             cdyn = max(cdyn,0.)   ! AS ALTERNATIVE, could clip ediff
>             ediff(i,1) = param(2)+cdyn*dg2(i,1)*snrm(i,1)
>          enddo
>       endif
>
> c     if (e.eq.nelv) call outpost(num,den,snrm,den,ediff,'dif')
> c     if (e.eq.nelv) call exitt
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine set_ds_filt(fh,fht,nt,diag,nx) ! setup test filter
>
>       INCLUDE 'SIZE'
>
>       real fh(nx*nx),fht(nx*nx),diag(nx)
>
> c Construct transfer function
>       call rone(diag,nx)
>
> c      diag(nx-0) = 0.01
> c      diag(nx-1) = 0.10
> c      diag(nx-2) = 0.50
> c      diag(nx-3) = 0.90
> c      diag(nx-4) = 0.99
> c      nt = nx - 2
>
>       diag(nx-0) = 0.05
>       diag(nx-1) = 0.50
>       diag(nx-2) = 0.95
>       nt = nx - 1
>
>       call build_1d_filt(fh,fht,diag,nx,nid)
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine planar_average_r(ua,u,w1,w2)
> c
> c     Compute s-t planar average of quantity u()
> c
>       include 'SIZE'
>       include 'GEOM'
>       include 'PARALLEL'
>       include 'WZ'
>       include 'ZPER'
> c
>       real ua(nx1,lelx),u(nx1,ny1,nx1,nelv),w1(nx1,lelx),w2(nx1,lelx)
>       integer e,eg,ex,ey,ez
> c
>       nx = nx1*nelx
>       call rzero(ua,nx)
>       call rzero(w1,nx)
> c
>       do e=1,nelt
> c
>          eg = lglel(e)
>          call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
> c
>          do i=1,nx1
>          do k=1,nz1
>          do j=1,ny1
>             zz = (1.-zgm1(i,1))/2.  ! = 1 for i=1, = 0 for k=nx1
>             aa = zz*area(j,k,4,e) + (1-zz)*area(j,k,2,e)  ! wgtd jacobian
>             w1(i,ex) = w1(i,ex) + aa
>             ua(i,ex) = ua(i,ex) + aa*u(i,j,k,e)
>          enddo
>          enddo
>          enddo
>       enddo
> c
>       call gop(ua,w2,'+  ',nx)
>       call gop(w1,w2,'+  ',nx)
> c
>       do i=1,nx
>          ua(i,1) = ua(i,1) / w1(i,1)   ! Normalize
>       enddo
> c
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine planar_average_s(ua,u,w1,w2)
> c
> c     Compute r-t planar average of quantity u()
> c
>       include 'SIZE'
>       include 'GEOM'
>       include 'PARALLEL'
>       include 'WZ'
>       include 'ZPER'
> c
>       real ua(ny1,nely),u(nx1,ny1,nx1,nelv),w1(ny1,nely),w2(ny1,nely)
>       integer e,eg,ex,ey,ez
> c
>       ny = ny1*nely
>       call rzero(ua,ny)
>       call rzero(w1,ny)
> c
>       do e=1,nelt
>          eg = lglel(e)
>          call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
> c
>          do k=1,nz1
>          do j=1,ny1
>          do i=1,nx1
>             zz = (1.-zgm1(j,2))/2.  ! = 1 for i=1, = 0 for k=nx1
>             aa = zz*area(i,k,1,e) + (1-zz)*area(i,k,3,e)  ! wgtd jacobian
>             w1(j,ey) = w1(j,ey) + aa
>             ua(j,ey) = ua(j,ey) + aa*u(i,j,k,e)
>          enddo
>          enddo
>          enddo
>       enddo
> c
>       call gop(ua,w2,'+  ',ny)
>       call gop(w1,w2,'+  ',ny)
> c
>       do i=1,ny
>          ua(i,1) = ua(i,1) / w1(i,1)   ! Normalize
>       enddo
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine planar_fill_s(u,ua)
> c
> c     Fill array u with planar values from ua().
> c     For tensor-product array of spectral elements
> c
>       include 'SIZE'
>       include 'GEOM'
>       include 'PARALLEL'
>       include 'WZ'
>       include 'ZPER'
>
>
>       real u(nx1,ny1,nz1,nelv),ua(ly1,lely)
>
>       integer e,eg,ex,ey,ez
>
>       melxyz = nelx*nely*nelz
>       if (melxyz.ne.nelgt) then
>          write(6,*) nid,' Error in planar_fill_s'
>      $                 ,nelgt,melxyz,nelx,nely,nelz
>          call exitt
>       endif
>
>       do e=1,nelt
>          eg = lglel(e)
>          call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
>
>          do j=1,ny1
>          do k=1,nz1
>          do i=1,nx1
>             u(i,j,k,e) = ua(j,ey)
>          enddo
>          enddo
>          enddo
>
>       enddo
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine set_grid_spacing(dg2)
> c
> c     Compute D^2, the grid spacing used in the DS sgs model.
> c
>       include 'SIZE'
>       include 'TOTAL'
>
>
>       real dg2(nx1,ny1,nz1,nelv)
>
>       integer e,eg,ex,ey,ez
>
>       gamma = 1.
>       gamma = gamma/ndim
>
>       n = nx1*ny1*nz1*nelv
>       call rone(dg2,n)
>       return               ! Comment this line for a non-trivial Delta defn
>
>       do e=1,nelv
>
>          do k=1,nz1
>            km = max(1  ,k-1)
>            kp = min(nz1,k+1)
>
>            do j=1,ny1
>              jm = max(1  ,j-1)
>              jp = min(ny1,j+1)
>
>              do i=1,nx1
>                im = max(1  ,i-1)
>                ip = min(nx1,i+1)
>
>                di = (xm1(ip,j,k,e)-xm1(im,j,k,e))**2
>      $            + (ym1(ip,j,k,e)-ym1(im,j,k,e))**2
>      $            + (zm1(ip,j,k,e)-zm1(im,j,k,e))**2
>
>                dj = (xm1(i,jp,k,e)-xm1(i,jm,k,e))**2
>      $            + (ym1(i,jp,k,e)-ym1(i,jm,k,e))**2
>      $            + (zm1(i,jp,k,e)-zm1(i,jm,k,e))**2
>
>                dk = (xm1(i,j,kp,e)-xm1(i,j,km,e))**2
>      $            + (ym1(i,j,kp,e)-ym1(i,j,km,e))**2
>      $            + (zm1(i,j,kp,e)-zm1(i,j,km,e))**2
>
>                di = di/(ip-im)
>                dj = dj/(jp-jm)
>                dk = dk/(kp-km)
>                dg2(i,j,k,e) = (di*dj*dk)**gamma
>
>              enddo
>            enddo
>          enddo
>       enddo
>
>       call dsavg(dg2)  ! average neighboring elements
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine wall_normal_average_s(u,ny,nel,v,w)
>       real u(ny,nel),w(1),v(1)
>       integer e
>
>       k=0
>       do e=1,nel    ! get rid of duplicated (ny,e),(1,e+1) points
>       do i=1,ny-1
>          k=k+1
>          w(k) = u(i,e)
>       enddo
>       enddo
>       k=k+1
>       w(k) = u(ny,nel)
>       n=k
>
>       npass = 2     ! Smooth
>       alpha = 0.2
>
>       do ipass=1,npass
>          do k=2,n-1
>             v(k) = (1.-alpha)*w(k) + 0.5*alpha*(w(k-1)+w(k+1))
>          enddo
>
>          do k=1,n
>             w(k) = v(k)
>          enddo
>       enddo
>
>       k=0
>       do e=1,nel    ! restore duplicated (ny,e),(1,e+1) points
>          do i=1,ny-1
>             k=k+1
>             u(i,e) = w(k)
>          enddo
>       enddo
>       k=k+1
>       u(ny,nel)=w(k)
>
>       do e=1,nel-1    ! restore duplicated (ny,e),(1,e+1) points
>          u(ny,e) = u(1,e+1)
>       enddo
>       return
>       end
> c-----------------------------------------------------------------------
> c  subroutines that follow are for fintpts based method
> c-----------------------------------------------------------------------
>       subroutine field_copy_si(fieldout,fieldin,idlist,nptsi)
>       include 'SIZE'
>       include 'TOTAL'
>
>       real    fieldin(1),fieldout(1)
>       integer idlist(1)
>
>       do i=1,nptsi
>         idx = idlist(i)
>         fieldout(idx) = fieldin(i)
>       enddo
>
>       return
>       end
> C--------------------------------------------------------------------------
>       subroutine field_eval_si(fieldout,fieldstride,fieldin)
>       include 'SIZE'
>       include 'TOTAL'
>
>       real fieldout(1),fieldin(1)
>
>       integer fieldstride,nptsi
>
>       parameter (lt=lelv*lx1*lz1)
>
>       integer elid_si(lt),proc_si(lt),ptid(lt),rcode_si(lt)
>       common /ptlist_int/ elid_si,proc_si,ptid,rcode_si,nptsi
>
>       real    rst_si(lt*ldim)
>       common /ptlist_real/ rst_si
>
>       integer inth_si
>       common / fpt_h_si/ inth_si
>
> c     Used for fgslib_findpts_eval of various fields
>       call fgslib_findpts_eval(inth_si,fieldout,fieldstride,
>      &                     rcode_si,1,
>      &                     proc_si,1,
>      &                     elid_si,1,
>      &                     rst_si,ndim,nptsi,
>      &                     fieldin)
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine rescale_inflow_fpt(ubar_in)  ! rescale inflow
>       include 'SIZE'
>       include 'TOTAL'
>
>       integer icalld,e,eg,f
>       save    icalld
>       data    icalld /0/
>       common /cvelbc/ uin(lx1,ly1,lz1,lelv)
>      $              , vin(lx1,ly1,lz1,lelv)
>      $              , win(lx1,ly1,lz1,lelv)
>
>       call get_flux_and_area(ubar,abar)
>       ubar  = ubar/abar        ! Ubar
>       scale = ubar_in/ubar     ! Scale factor
>
>       if (nid.eq.0.and.(istep.le.100.or.mod(istep,100).eq.0))
>      $  write(6,1) istep,time,scale,ubar,abar
>     1   format(1i8,1p4e14.6,' rescale')
>
> c     Rescale the flow to match ubar_in
>       do e=1,nelv
>       do f=1,2*ldim
>         if (cbc(f,e,1).eq.'v  ') then
>            call facind (kx1,kx2,ky1,ky2,kz1,kz2,nx1,ny1,nz1,f)
>            do iz=kz1,kz2
>            do iy=ky1,ky2
>            do ix=kx1,kx2
>               uin(ix,iy,iz,e) =  scale*uin(ix,iy,iz,e)
>               vin(ix,iy,iz,e) =  scale*vin(ix,iy,iz,e)
>               win(ix,iy,iz,e) =  scale*win(ix,iy,iz,e)
>            enddo
>            enddo
>            enddo
>         endif
>       enddo
>       enddo
>
>       ifield = 1       ! Project into H1, just to be sure....
>       call dsavg(uin)
>       call dsavg(vin)
>       if (ldim.eq.3) call dsavg(win)
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine get_flux_and_area(vvflux,vvarea)
>       include 'SIZE'
>       include 'TOTAL'
>       common /cvelbc/ uin(lx1,ly1,lz1,lelv)
>      $              , vin(lx1,ly1,lz1,lelv)
>      $              , win(lx1,ly1,lz1,lelv)
>       real vvflux,vvarea
>       real work(lx1*ly1*lz1)
>       integer e,f
>
>       nxz   = nx1*nz1
>       nface = 2*ndim
>
>       vvflux = 0.
>       vvarea = 0.
>
>       do e=1,nelv
>       do f=1,nface
>          if (cbc(f,e,1).eq.'v  ') then
>             call surface_flux(dq,uin,vin,win,e,f,work)
>             vvflux = vvflux + dq
>             vvarea = vvarea + vlsum(area(1,1,f,e),nxz)
>          endif
>       enddo
>       enddo
>       vvflux = glsum(vvflux,1)
>       vvarea = glsum(vvarea,1)
>       vvflux = -vvflux !flux in is negative
>
>       return
>       end
> c-----------------------------------------------------------------------
>       subroutine set_inflow_fpt_setup(dxx,dyy,dzz)   ! set up inflow BCs
>       include 'SIZE'
>       include 'TOTAL'
> c
> c setup recirculation boundary condition based on user supplied dx,dy,dz
> c dx,dy,dz is the vector from the inflow where the user wants the velocity
> c data to be interpolated from
> c
>       integer icalld,e,eg,i,f,nptsi
>       save    icalld
>       data    icalld /0/
>       real dxx,dyy,dzz
>
>       parameter (lt=lx1*lz1*lelv)
>       real rst_si(lt*ldim),xyz_si(lt*ldim)
>       real dist_si(lt),vals_si(lt)
>
>       integer elid_si(lt), proc_si(lt),ptid(lt)
>       integer rcode_si(lt)
>       common /ptlist_real/ rst_si
>       common /ptlist_int/ elid_si,proc_si,ptid,rcode_si,nptsi
>       integer inth_si
>       common / fpt_h_si/ inth_si
>       common /cvelbc/ uin(lx1,ly1,lz1,lelv)
>      $              , vin(lx1,ly1,lz1,lelv)
>      $              , win(lx1,ly1,lz1,lelv)
>       common /nekmpi/ nidd,npp,nekcomm,nekgroup,nekreal
>
>       n = nx1*ny1*nz1*nelv
> ccc
> c     Gather info for findpts
> ccc
>       nptsi = 0
>       nxyz = nx1*ny1*nz1
>
>       do e=1,nelv
>       do f=1,2*ndim  !Identify the xyz of the points that are to be found
>        if (cbc(f,e,1).eq.'v  ') then
>            call facind (kx1,kx2,ky1,ky2,kz1,kz2,nx1,ny1,nz1,f)
>            do iz=kz1,kz2
>            do iy=ky1,ky2
>            do ix=kx1,kx2
>             nptsi = nptsi+1
>             xyz_si(ldim*(nptsi-1)+1) = xm1(ix,iy,iz,e) + dxx
>             xyz_si(ldim*(nptsi-1)+2) = ym1(ix,iy,iz,e) + dyy
>       if (ldim.eq.3) xyz_si(ldim*(nptsi-1)+ldim) = zm1(ix,iy,iz,e) + dzz
>             ptid(nptsi) = (e-1)*nxyz+(iz-1)*lx1*ly1+(iy-1)*lx1+ix
>            enddo
>            enddo
>            enddo
>        endif
>       enddo
>       enddo
>       mptsi=iglmax(nptsi,1)
>       if (mptsi.gt.lt)
>      $  call exitti('ERROR: increase lt in inflow_fpt routines.$',mptsi)
>
> c     Setup findpts
>
>       tol     = 1e-10
>       npt_max = 256
>       nxf     = 2*nx1 ! fine mesh for bb-test
>       nyf     = 2*ny1
>       nzf     = 2*nz1
>       bb_t    = 0.1 ! relative size to expand bounding boxes by
>       bb_t    = 0.1 ! relative size to expand bounding boxes by
>       call fgslib_findpts_setup(inth_si,nekcomm,npp,ndim,
>      &                   xm1,ym1,zm1,nx1,ny1,nz1,
>      &                   nelt,nxf,nyf,nzf,bb_t,n,n,
>      &                   npt_max,tol)
>
>
> c     Call findpts to determine el,proc,rst of the xyzs determined above
>
>       call fgslib_findpts(inth_si,rcode_si,1,
>      &             proc_si,1,
>      &             elid_si,1,
>      &             rst_si,ndim,
>      &             dist_si,1,
>      &             xyz_si(1),ldim,
>      &             xyz_si(2),ldim,
>      &             xyz_si(3),ldim,nptsi)
>
>       return
>       end
> C-----------------------------------------------------------------------
>       subroutine set_inflow_fpt(dxx,dyy,dzz,ubar)   ! set up inflow BCs
>       include 'SIZE'
>       include 'TOTAL'
>
> c setup recirculation boundary condition based on user supplied dx,dy,dz
> c dx,dy,dz is the vector from the inflow where the user wants the
> c velocity data to be interpolated from
>
>       integer icalld
>       save    icalld
>       data    icalld /0/
>       real dxx,dyy,dzz
>
>       parameter (lt=lx1*lz1*lelv)
>       real rst_si(lt*ldim),xyz_si(lt*ldim)
>       real dist_si(lt),vals_si(lt)
>       common /ptlist_real/ rst_si
>
>       integer elid_si(lt), proc_si(lt),ptid(lt),rcode_si(lt)
>       common /ptlist_int/ elid_si,proc_si,ptid,rcode_si,nptsi
>       integer inth_si
>       common / fpt_h_si/ inth_si
>       common /cvelbc/ uin(lx1,ly1,lz1,lelv)
>      $              , vin(lx1,ly1,lz1,lelv)
>      $              , win(lx1,ly1,lz1,lelv)
>
>
> c     Gather info for findpts and set up inflow BC
>       if (icalld.eq.0) call set_inflow_fpt_setup(dxx,dyy,dzz)
>       icalld=1
>
>
> c     Eval fields and copy to uvwin array
>       call field_eval_si(vals_si,1,vx)
>       call field_copy_si(uin,vals_si,ptid,nptsi)
>
>       call field_eval_si(vals_si,1,vy)
>       call field_copy_si(vin,vals_si,ptid,nptsi)
>
>       if (ldim.eq.3) then
>         call field_eval_si(vals_si,1,vz)
>         call field_copy_si(win,vals_si,ptid,nptsi)
>       endif
>
> c     Rescale the flow so that ubar,vbar or wbar is ubar
>       call rescale_inflow_fpt(ubar)
>
>       return
>       end
> C-----------------------------------------------------------------------
>
> c automatically added by makenek
>       subroutine usrsetvert(glo_num,nel,nx,ny,nz) ! to modify glo_num
>       integer*8 glo_num(1)
>       return
>       end
>
> ------------------------------
>
> Subject: Digest Footer
>
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>
>
> ------------------------------
>
> End of Nek5000-users Digest, Vol 112, Issue 10
> **********************************************
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20180607/467f9679/attachment-0001.html>


More information about the Nek5000-users mailing list