<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>