[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