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