[Nek5000-users] Extruding 2D .fld file to 3D
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Thu Oct 24 06:40:16 CDT 2013
Hi Paul,
This works perfectly! Thanks for your help!
Holly JOHNSON
PhD Student ONERA - DAFE
Département d'Aérodynamique Fondamentale et Expérimentale
8, Rue des Vertugadins
92190 - Meudon - France
Le 24.10.2013 02:52, nek5000-users at lists.mcs.anl.gov a écrit :
> Hi Holly,
>
>
> I would likely do this in usrdat2().
>
> Suppose my 2D problem had 70 elements. The the following
> should work:
>
> In SIZE, set lelx=70,lely=1, and lelz=# of levels in your extruded
> geometry.
>
> In usrdat2() in the .usr file, set nelx=lelx,nely=lely,nelz=lelz.
>
> Then:
> parameter (l2d=lx1*ly1*lelx) ! Number of points in 2D field
> common /mystuff/ u2d(l2d),v2d(l2d)
>
> if (nid.eq.0) then
> open(33,file='my2dfile.dat')
> do i=1,l2d
> read(33,*) u2d(i),v2d(i)
> enddo
> close(33)
> endif
> nbytes = 8*l2d
> call bcast(u2d,nbytes)
> call bcast(v2d,nbytes)
>
> Once you load your 2D data you can use
>
> call z_average_transpose(vx,u2d)
> call z_average_transpose(vy,v2d)
>
> to map this to 3D.
>
> After that, I would just dump the 3D results and exit:
>
> call outpost(vx,vy,vz,pr,t,' ')
> call exitti('Quit in usrdat2.$',nelgt)
>
> This _should_ work, but probably needs a couple iterations
> to get everything correct.
>
> To generate myfile2d.dat, I would do the following:
>
> 1) Run the 2D code and dump an ascii .fld file by setting
> param(66)=0 in usrdat2 or userchk prior to output the fld file.
>
> 2) Edit the resultant .fld file and get rid of the first few
> line at the top of that file so that your u,v data is all
> that it holds. (Look both at the top and bottom of the file
> when you're ready -- they should have the same column format.)
>
>
> Hope this helps.
>
> Paul
>
>
>
>
>
> On Wed, 23 Oct 2013, nek5000-users at lists.mcs.anl.gov wrote:
>
>> Hello Neks !
>>
>> In order to reduce runtime, it would be useful to be able to run a 2D
>> simulation and extrude the results to a 3D domain when using a 2D
>> model (for example the Rankine vortex). Is it possible to extrude a 2D
>> .fld file to a 3D one ? I was thinking of loading a 2D .fld file into
>> a 3D simulation using something similar to load_fld() and then using
>> loops to impose the same velocity to all the points in the domain.
>> However load_fld() doesn't seem to like loading a 2D .fld in a 3D
>> simulation. Has anyone else had this problem ?
>>
>> Thanks for your help!
>>
>> Holly JOHNSON
>>
>> PhD Student ONERA - DAFE
>> Département d'Aérodynamique Fondamentale et Expérimentale
>> 8, Rue des Vertugadins
>> 92190 - Meudon - France
>>
>> _______________________________________________
>> Nek5000-users mailing list
>> Nek5000-users at lists.mcs.anl.gov
>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>>
> uuuu
> c-----------------------------------------------------------------------
> subroutine q_filter(wght)
> c
> c filter vx,vy,vz, and p by simple interpolation
> c
> include 'SIZE'
> include 'TOTAL'
> c
> c
> c These are the dimensions that we interpolate onto for v and p:
> parameter(lxv=lx1-1)
> parameter(lxp=lx2-1)
> c
> real intdv(lx1,lx1)
> real intuv(lx1,lx1)
> real intdp(lx1,lx1)
> real intup(lx1,lx1)
> real intv(lx1,lx1)
> real intp(lx1,lx1)
> c
> save intdv
> save intuv
> save intdp
> save intup
> save intv
> save intp
>
> common /ctmp0/ intw,intt
> $ , wk1,wk2
> $ , zgmv,wgtv,zgmp,wgtp,tmax(100),omax(103)
>
> real intw(lx1,lx1)
> real intt(lx1,lx1)
> real wk1 (lx1,lx1,lx1,lelt)
> real wk2 (lx1,lx1,lx1)
> real zgmv(lx1),wgtv(lx1),zgmp(lx1),wgtp(lx1)
> c
> c outpost arrays
> parameter (lt=lx1*ly1*lz1*lelv)
> common /scruz/ w1(lt),w2(lt),w3(lt),wt(lt)
>
> character*18 sfmt
>
> integer icalld
> save icalld
> data icalld /0/
>
>
> imax = nid
> imax = iglmax(imax,1)
> jmax = iglmax(imax,1)
> if (icalld.eq.0) then
> icalld = 1
> ncut = param(101)+1
> call build_new_filter(intv,zgm1,nx1,ncut,wght,nid)
> elseif (icalld.lt.0) then ! old (std.) filter
> icalld = 1
> call zwgll(zgmv,wgtv,lxv)
> call igllm(intuv,intw,zgmv,zgm1,lxv,nx1,lxv,nx1)
> call igllm(intdv,intw,zgm1,zgmv,nx1,lxv,nx1,lxv)
> c
> call zwgl (zgmp,wgtp,lxp)
> call iglm (intup,intw,zgmp,zgm2,lxp,nx2,lxp,nx2)
> call iglm (intdp,intw,zgm2,zgmp,nx2,lxp,nx2,lxp)
> c
> c Multiply up and down interpolation into single operator
> c
> call mxm(intup,nx2,intdp,lxp,intp,nx2)
> call mxm(intuv,nx1,intdv,lxv,intv,nx1)
> c
> c Weight the filter to make it a smooth (as opposed to
> truncated)
> c decay in wave space
> c
> w0 = 1.-wght
> call ident(intup,nx2)
> call add2sxy(intp,wght,intup,w0,nx2*nx2)
> c
> call ident (intuv,nx1)
> call add2sxy (intv ,wght,intuv,w0,nx1*nx1)
> c
> c if (nid.eq.0) call outmatx(intp,nx2,nx2,21,'flt2')
> c if (nid.eq.0) call outmatx(zgm2 ,nx2,1 ,22,'zgm2')
> c if (nid.eq.0) call outmatx(intv,nx1,nx1,11,'flt1')
> c if (nid.eq.0) call outmatx(zgm1 ,nx1,1 ,12,'zgm1')
> c
> endif
>
> ifldt = ifield
> c ifield = 1
>
> if ( (ifflow.and. .not. ifmhd) .or.
> $ (ifield.eq.1 .and. ifmhd) ) then
> call filterq(vx,intv,nx1,nz1,wk1,wk2,intt,if3d,umax)
> call filterq(vy,intv,nx1,nz1,wk1,wk2,intt,if3d,vmax)
> if (if3d)
> $ call filterq(vz,intv,nx1,nz1,wk1,wk2,intt,if3d,wmax)
> if (ifsplit.and..not.iflomach)
> $ call filterq(pr,intv,nx1,nz1,wk1,wk2,intt,if3d,pmax)
> endif
> c
> if (ifmhd.and.ifield.eq.ifldmhd) then
> call filterq(bx,intv,nx1,nz1,wk1,wk2,intt,if3d,umax)
> call filterq(by,intv,nx1,nz1,wk1,wk2,intt,if3d,vmax)
> if (if3d)
> $ call filterq(bz,intv,nx1,nz1,wk1,wk2,intt,if3d,wmax)
> endif
> c
> if (ifpert) then
> do j=1,npert
>
> ifield = 1
> call filterq(vxp(1,j),intv,nx1,nz1,wk1,wk2,intt,if3d,umax)
> call filterq(vyp(1,j),intv,nx1,nz1,wk1,wk2,intt,if3d,vmax)
> if (if3d)
> $ call filterq(vzp(1,j),intv,nx1,nz1,wk1,wk2,intt,if3d,wmax)
>
> ifield = 2
> if (ifheat .and. .not.ifcvode)
> $ call filterq(tp(1,j,1),intv,nx1,nz1,wk1,wk2,intt,if3d,wmax)
>
> enddo
> endif
> c
> mmax = 0
> if (ifflow) then
> c pmax = glmax(pmax,1)
> omax(1) = glmax(umax,1)
> omax(2) = glmax(vmax,1)
> omax(3) = glmax(wmax,1)
> mmax = ndim
> endif
>
> c
> nfldt = 1+npscal
> if (ifheat .and. .not.ifcvode) then
> do ifld=1,nfldt
> ifield = ifld + 1
> call filterq(t(1,1,1,1,ifld),intv
> $ ,nx1,nz1,wk1,wk2,intt,if3d,tmax(ifld))
> mmax = mmax+1
> omax(mmax) = glmax(tmax(ifld),1)
> enddo
> endif
>
> if (nid.eq.0) then
> if (npscal.eq.0) then
> c write(6,101) mmax
> c write(sfmt,101) mmax
> c 101 format('''(i8,1p',i1,'e12.4,a6)''')
> c write(6,sfmt) istep,(omax(k),k=1,mmax),' qfilt'
> c write(6,'(i8,1p4e12.4,a6)') istep,(omax(k),k=1,mmax),' qfilt'
> else
> if (if3d) then
> write(6,1) istep,ifield,umax,vmax,wmax
> else
> write(6,1) istep,ifield,umax,vmax
> endif
> 1 format(4x,i7,i3,' qfilt:',1p3e12.4)
> if(ifheat .and. .not.ifcvode)
> & write(6,'(1p50e12.4)') (tmax(k),k=1,nfldt)
> endif
> endif
>
> ifield = ifldt ! RESTORE ifield
>
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine filterq(v,f,nx,nz,w1,w2,ft,if3d,dmax)
> c
> include 'SIZE'
> include 'TSTEP'
>
> real v(nx*nx*nz,nelt),w1(1),w2(1)
> logical if3d
> c
> real f(nx,nx),ft(nx,nx)
> c
> integer e
> c
> call transpose(ft,nx,f,nx)
> c
> nxyz=nx*nx*nz
> dmax = 0.
>
>
> nel = nelfld(ifield)
>
>
> if (if3d) then
> do e=1,nel
> c Filter
> call copy(w2,v(1,e),nxyz)
> call mxm(f,nx,w2,nx,w1,nx*nx)
> i=1
> j=1
> do k=1,nx
> call mxm(w1(i),nx,ft,nx,w2(j),nx)
> i = i+nx*nx
> j = j+nx*nx
> enddo
> call mxm (w2,nx*nx,ft,nx,w1,nx)
> call sub3(w2,v(1,e),w1,nxyz)
> call copy(v(1,e),w1,nxyz)
> smax = vlamax(w2,nxyz)
> dmax = max(dmax,abs(smax))
> enddo
> c
> else
> do e=1,nel
> c Filter
> call copy(w1,v(1,e),nxyz)
> call mxm(f ,nx,w1,nx,w2,nx)
> call mxm(w2,nx,ft,nx,w1,nx)
> c
> call sub3(w2,v(1,e),w1,nxyz)
> call copy(v(1,e),w1,nxyz)
> smax = vlamax(w2,nxyz)
> dmax = max(dmax,abs(smax))
> enddo
> endif
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine outmatx(a,m,n,io,name)
> real a(m*n)
> character*4 name
> c
> open(unit=io,file=name)
> do i=1,m*n
> write(io,1) a(i)
> enddo
> 1 format(1p1e22.13)
> close(unit=io)
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine drag_calc(scale)
> c
> INCLUDE 'SIZE'
> INCLUDE 'TOTAL' c
> common /scrns/ pm1(lx1,ly1,lz1,lelv)
>
> $,vxx(lx1,ly1,lz1,lelv),vxy(lx1,ly1,lz1,lelv),vxz(lx1,ly1,lz1,lelv)
>
> $,vyx(lx1,ly1,lz1,lelv),vyy(lx1,ly1,lz1,lelv),vyz(lx1,ly1,lz1,lelv)
> common /scruz/
> $
> vzx(lx1,ly1,lz1,lelv),vzy(lx1,ly1,lz1,lelv),vzz(lx1,ly1,lz1,lelv)
> $,one(lx1,ly1,lz1,lelv)
> real work(1)
> equivalence (work,one)
> c
> common /cdrag/ dragx(0:maxobj),dragy(0:maxobj),dragz(0:maxobj)
> $ , momx(0:maxobj), momy(0:maxobj), momz(0:maxobj)
> $ , dpdx_mean,dpdy_mean,dpdz_mean
> real momx ,momy ,momz
> c
> common /tdrag/ drag(11)
> real dragpx,dragpy,dragpz,dragvx,dragvy,dragvz
> real momvx ,momvy ,momvz
> real check1,check2
> c
> equivalence (dragpx,drag(1)),(dragpy,drag(2)),(dragpz,drag(3))
> equivalence (dragvx,drag(4)),(dragvy,drag(5)),(dragvz,drag(6))
> equivalence (momvx ,drag(7)),(momvy ,drag(8)),(momvz ,drag(9))
> equivalence (check1,drag(10)),(check2,drag(11))
>
> common /cvflow_r/ flow_rate,base_flow,domain_length,xsec
> $ , scale_vf(3)
>
> ntot1 = nx1*ny1*nz1*nelv
>
> c Map pressure onto mesh 1 (vxx and vyy are used as work arrays)
> call mappr(pm1,pr,vxx,vyy)
> call rone (one,ntot1)
> c
> c Add mean_pressure_gradient.X to p:
>
> if (param(55).ne.0) then
> dpdx_mean = -scale_vf(1)
> dpdy_mean = -scale_vf(2)
> dpdz_mean = -scale_vf(3)
> endif
>
> call add2s2(pm1,xm1,dpdx_mean,ntot1) ! Doesn't work if object is
> cut by
> call add2s2(pm1,ym1,dpdy_mean,ntot1) ! periodicboundary. In
> this case,
> call add2s2(pm1,zm1,dpdz_mean,ntot1) ! set ._mean=0 and
> compensate in
> c ! usrchk() [ pff 10/21/04
> ].
>
> c Compute du/dn
> CALL DUDXYZ (vxx,vx,RXM1,SXM1,TXM1,JACM1,1,1)
> CALL DUDXYZ (vxy,vx,RYM1,SYM1,TYM1,JACM1,1,1)
> if (if3d) CALL DUDXYZ (vxz,vx,RZM1,SZM1,TZM1,JACM1,1,1)
> c
> CALL DUDXYZ (vyx,vy,RXM1,SXM1,TXM1,JACM1,1,1)
> CALL DUDXYZ (vyy,vy,RYM1,SYM1,TYM1,JACM1,1,1)
> if (if3d) CALL DUDXYZ (vyz,vy,RZM1,SZM1,TZM1,JACM1,1,1)
> c
> if (if3d) then
> CALL DUDXYZ (vzx,vz,RXM1,SXM1,TXM1,JACM1,1,1)
> CALL DUDXYZ (vzy,vz,RYM1,SYM1,TYM1,JACM1,1,1)
> CALL DUDXYZ (vzz,vz,RZM1,SZM1,TZM1,JACM1,1,1)
> endif
> c
> c Fill up viscous array w/ default
> c
> if (istep.lt.1) call cfill(vdiff,param(2),ntot1)
> c
> call col2(vxx,vdiff,ntot1)
> call col2(vxy,vdiff,ntot1)
> call col2(vxz,vdiff,ntot1)
> call col2(vyx,vdiff,ntot1)
> call col2(vyy,vdiff,ntot1)
> call col2(vyz,vdiff,ntot1)
> call col2(vzx,vdiff,ntot1)
> call col2(vzy,vdiff,ntot1)
> call col2(vzz,vdiff,ntot1)
> c
> dragxt=0.0
> dragyt=0.0
> dragzt=0.0
> c
> DO 100 II=1,NHIS
> IF (HCODE(10,II).NE.'I') GOTO 100
> IOBJ = LOCHIS(1,II)
> MEMTOT = NMEMBER(IOBJ)
> C
> c
> IF (HCODE(1,II).NE.' ' .OR. HCODE(2,II).NE.' ' .OR.
> $ HCODE(3,II).NE.' ' ) THEN
> IFIELD = 1
> c
> c---------------------------------------------------------------------------
> c Compute drag for this object
> c---------------------------------------------------------------------------
> c
> dragvx=0.0
> dragvy=0.0
> dragvz=0.0
> dragpx=0.0
> dragpy=0.0
> dragpz=0.0
> c
> momvx=0.0
> momvy=0.0
> momvz=0.0
> c
> check1=0.0
> check2=0.0
> DO 50 MEM=1,MEMTOT
> ISK = 0
> IEG = OBJECT(IOBJ,MEM,1)
> IFC = OBJECT(IOBJ,MEM,2)
> IF (GLLNID(IEG).EQ.NID) THEN
> C This processor has a contribution
> IE = GLLEL(IEG)
> c
> c Pressure drag
> check1=check1+facint(one,one,area,ifc,ie)
> check2=check2+facint(one,uny,area,ifc,ie)
> c
> dragpx=dragpx+facint(pm1,unx,area,ifc,ie)
> dragpy=dragpy+facint(pm1,uny,area,ifc,ie)
> if (if3d) dragpz=dragpz+facint(pm1,unz,area,ifc,ie)
> c
> c Viscous drag
> if (if3d) then
> dragvx=dragvx+facint(vxx,unx,area,ifc,ie)
> $ +facint(vxy,uny,area,ifc,ie)
> $ +facint(vxz,unz,area,ifc,ie)
> $ +facint(vxx,unx,area,ifc,ie)
> $ +facint(vyx,uny,area,ifc,ie)
> $ +facint(vzx,unz,area,ifc,ie)
> dragvy=dragvy+facint(vyx,unx,area,ifc,ie)
> $ +facint(vyy,uny,area,ifc,ie)
> $ +facint(vyz,unz,area,ifc,ie)
> $ +facint(vxy,unx,area,ifc,ie)
> $ +facint(vyy,uny,area,ifc,ie)
> $ +facint(vzy,unz,area,ifc,ie)
> dragvz=dragvz+facint(vzx,unx,area,ifc,ie)
> $ +facint(vzy,uny,area,ifc,ie)
> $ +facint(vzz,unz,area,ifc,ie)
> $ +facint(vxz,unx,area,ifc,ie)
> $ +facint(vyz,uny,area,ifc,ie)
> $ +facint(vzz,unz,area,ifc,ie)
> c
> momvx=momvx-facint2(vxy,unx,unz,area,ifc,ie)
> $ -facint2(vyx,unx,unz,area,ifc,ie)
> $ -facint2(vyy,uny,unz,area,ifc,ie)
> $ -facint2(vyy,uny,unz,area,ifc,ie)
> $ -facint2(vzy,unz,unz,area,ifc,ie)
> $ -facint2(vyz,unz,unz,area,ifc,ie)
> $ +facint2(vxz,unx,uny,area,ifc,ie)
> $ +facint2(vzx,unx,uny,area,ifc,ie)
> $ +facint2(vyz,uny,uny,area,ifc,ie)
> $ +facint2(vzy,uny,uny,area,ifc,ie)
> $ +facint2(vzz,unz,uny,area,ifc,ie)
> $ +facint2(vzz,unz,uny,area,ifc,ie)
> momvy=momvy+facint2(vxx,unx,unz,area,ifc,ie)
> $ +facint2(vxx,unx,unz,area,ifc,ie)
> $ +facint2(vyx,uny,unz,area,ifc,ie)
> $ +facint2(vxy,uny,unz,area,ifc,ie)
> $ +facint2(vzx,unz,unz,area,ifc,ie)
> $ +facint2(vxz,unz,unz,area,ifc,ie)
> $ -facint2(vxz,unx,unx,area,ifc,ie)
> $ -facint2(vzx,unx,unx,area,ifc,ie)
> $ -facint2(vyz,uny,unx,area,ifc,ie)
> $ -facint2(vzy,uny,unx,area,ifc,ie)
> $ -facint2(vzz,unz,unx,area,ifc,ie)
> $ -facint2(vzz,unz,unx,area,ifc,ie)
> momvz=momvz-facint2(vxx,unx,uny,area,ifc,ie)
> $ -facint2(vxx,unx,uny,area,ifc,ie)
> $ -facint2(vyx,uny,uny,area,ifc,ie)
> $ -facint2(vxy,uny,uny,area,ifc,ie)
> $ -facint2(vzx,unz,uny,area,ifc,ie)
> $ -facint2(vxz,unz,uny,area,ifc,ie)
> $ +facint2(vxy,unx,unx,area,ifc,ie)
> $ +facint2(vyx,unx,unx,area,ifc,ie)
> $ +facint2(vyy,uny,unx,area,ifc,ie)
> $ +facint2(vyy,uny,unx,area,ifc,ie)
> $ +facint2(vzy,unz,unx,area,ifc,ie)
> $ +facint2(vyz,unz,unx,area,ifc,ie) c
> else
> dragvx=dragvx+facint(vxx,unx,area,ifc,ie)
> $ +facint(vxy,uny,area,ifc,ie)
> dragvy=dragvy+facint(vyx,unx,area,ifc,ie)
> $ +facint(vyy,uny,area,ifc,ie)
> endif
> ENDIF
> 50 CONTINUE
> c
> c Sum contributions from all processors
> call gop(drag,work,'+ ',11)
> dragvx = -dragvx
> dragvy = -dragvy
> dragvz = -dragvz
> ENDIF
> c
> c Scale by user specified scale factor (for convenience)
> c
> dragvx = scale*dragvx
> dragvy = scale*dragvy
> dragvz = scale*dragvz
> c
> dragpx = scale*dragpx
> dragpy = scale*dragpy
> dragpz = scale*dragpz
> c
> dragx(iobj) = dragvx+dragpx
> dragy(iobj) = dragvy+dragpy
> dragz(iobj) = dragvz+dragpz
> c
> c
> momx(iobj) = 0.5*momvx
> momy(iobj) = 0.5*momvy
> momz(iobj) = 0.5*momvz
> c
> dragxt = dragxt + dragx(iobj)
> dragyt = dragyt + dragy(iobj)
> dragzt = dragzt + dragz(iobj)
> c
> if (nid.eq.0.and.istep.eq.1)
> $ write(6,*) 'drag_calc: scale=',scale
> if (nid.eq.0) then
> write(6,6)
> istep,time,dragx(iobj),dragpx,dragvx,'dragx',iobj
> write(6,6)
> istep,time,dragy(iobj),dragpy,dragvy,'dragy',iobj
> if (if3d)
> $ write(6,6)
> istep,time,dragz(iobj),dragpz,dragvz,'dragz',iobj
> c
> c done by zly (3/17/03)
> c if(if3d) then
> c write(6,113) istep,time,momx,momy,momz
> c else
> c write(6,112) istep,time,momx,momy
> c endif
> c
> endif
> 6 format(i8,1p4e15.7,2x,a5,i5)
> 112 format(i6,1p3e15.7,' momx')
> 113 format(i6,1p4e15.7,' momx')
> if (istep.lt.10.and.nid.eq.0)
> $ write(6,9) 'check:',check1,check2,dpdx_mean,istep
> 9 format(a6,1p3e16.8,i9)
> c if (time.gt.1.0.and.dragx.gt.10.) call emerxit
> 100 continue
> c
> if (nid.eq.0) then
> write(6,6) istep,time,dragxt,dragpx,dragvx,'drgxt',iobj
> write(6,6) istep,time,dragyt,dragpy,dragvy,'drgyt',iobj
> if (if3d)
> $ write(6,6) istep,time,dragzt,dragpz,dragvz,'drgzt',iobj
> endif
> c
> dragx(0) = dragxt
> dragy(0) = dragyt
> dragz(0) = dragzt
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine mappr(pm1,pm2,pa,pb)
> c
> INCLUDE 'SIZE'
> INCLUDE 'TOTAL'
> real pm1(lx1,ly1,lz1,lelv),pm2(lx2,ly2,lz2,lelv)
> $ ,pa (lx1,ly2,lz2) ,pb (lx1,ly1,lz2)
> c
> C Map the pressure onto the velocity mesh
> C
> NGLOB1 = NX1*NY1*NZ1*NELV
> NYZ2 = NY2*NZ2
> NXY1 = NX1*NY1
> NXYZ = NX1*NY1*NZ1
> C
> IF (IFSPLIT) THEN
> CALL COPY(PM1,PM2,NGLOB1)
> ELSE
> DO 1000 IEL=1,NELV
> CALL MXM (IXM21,NX1,PM2(1,1,1,IEL),NX2,pa (1,1,1),NYZ2)
> DO 100 IZ=1,NZ2
> CALL MXM (PA(1,1,IZ),NX1,IYTM21,NY2,PB(1,1,IZ),NY1)
> 100 CONTINUE
> CALL MXM (PB(1,1,1),NXY1,IZTM21,NZ2,PM1(1,1,1,IEL),NZ1)
> 1000 CONTINUE
>
> C Average the pressure on elemental boundaries
> IFIELD=1
> CALL DSSUM (PM1,NX1,NY1,NZ1)
> CALL COL2 (PM1,VMULT,NGLOB1)
>
> ENDIF
> C
> C
> return
> end
> c
> c-----------------------------------------------------------------------
> function facint_a(a,area,f,e)
> c Integrate areal array a() on face f of element e. 27 June, 2012
> pff
>
> c f is in the preprocessor notation
>
> include 'SIZE'
> include 'TOPOL'
> real a(lx1,lz1,6,lelt),area(lx1,lz1,6,lelt)
>
> integer e,f
>
> sum=0.0
> do i=1,lx1*lz1
> sum = sum + a(i,1,f,e)*area(i,1,f,e)
> enddo
>
> facint_a = sum
>
> return
> end
> c-----------------------------------------------------------------------
> function facint_v(a,area,f,e)
> c Integrate volumetric array a() on face f of element e
>
> c f is in the preprocessor notation
> c fd is the dssum notation.
> c 27 June, 2012 PFF
>
> include 'SIZE'
> include 'TOPOL'
> real a(lx1,ly1,lz1,lelt),area(lx1,lz1,6,lelt)
>
> integer e,f,fd
>
> call dsset(nx1,ny1,nz1) ! set counters
> fd = eface1(f)
> js1 = skpdat(1,fd)
> jf1 = skpdat(2,fd)
> jskip1 = skpdat(3,fd)
> js2 = skpdat(4,fd)
> jf2 = skpdat(5,fd)
> jskip2 = skpdat(6,fd)
>
> sum=0.0
> i = 0
> do 100 j2=js2,jf2,jskip2
> do 100 j1=js1,jf1,jskip1
> i = i+1
> sum = sum + a(j1,j2,1,e)*area(i,1,f,e)
> 100 continue
>
> facint_v = sum
>
> return
> end
> c-----------------------------------------------------------------------
> function facint(a,b,area,ifc,ie)
> c
> C
> C Take the dot product of A and B on the surface IFACE1 of element
> IE.
> C
> C IFACE1 is in the preprocessor notation
> C IFACE is the dssum notation.
> C 5 Jan 1989 15:12:22 PFF
> C
> INCLUDE 'SIZE'
> INCLUDE 'TOPOL'
> DIMENSION A (LX1,LY1,LZ1,lelv)
> $ ,B (lx1,lz1,6,lelv)
> $ ,area (lx1,lz1,6,lelv)
> C
> C Set up counters
> C
> CALL DSSET(NX1,NY1,NZ1)
> IFACE = EFACE1(IFC)
> JS1 = SKPDAT(1,IFACE)
> JF1 = SKPDAT(2,IFACE)
> JSKIP1 = SKPDAT(3,IFACE)
> JS2 = SKPDAT(4,IFACE)
> JF2 = SKPDAT(5,IFACE)
> JSKIP2 = SKPDAT(6,IFACE)
> C
> SUM=0.0
> I = 0
> DO 100 J2=JS2,JF2,JSKIP2
> DO 100 J1=JS1,JF1,JSKIP1
> I = I+1
> SUM = SUM + A(J1,J2,1,ie)*B(I,1,ifc,ie)*area(I,1,ifc,ie)
> c SUM = SUM + A(J1,J2,1,ie)*B(J1,J2,1,ie)*area(I,1,ifc,ie)
> 100 CONTINUE
> C
> facint = SUM
> c
> return
> end
> c-----------------------------------------------------------------------
> function facint2(a,b,c,area,ifc,ie)
> include 'SIZE'
> include 'TOPOL'
> dimension a (lx1,ly1,lz1,lelv)
> $ , b (lx1,lz1,6,lelv)
> $ , c (lx1,lz1,6,lelv)
> $ , area (lx1,lz1,6,lelv)
> call dsset(nx1,ny1,nz1)
> iface = eface1(ifc)
> js1 = skpdat(1,iface)
> jf1 = skpdat(2,iface)
> jskip1 = skpdat(3,iface)
> js2 = skpdat(4,iface)
> jf2 = skpdat(5,iface)
> jskip2 = skpdat(6,iface)
> sum=0.0
> i=0
> do j2=js2,jf2,jskip2
> do j1=js1,jf1,jskip1
> i=i+1
> sum=sum+a(j1,j2,1,ie)*b(i,1,ifc,ie)*c(i,1,ifc,ie)
> $ *area(i,1,ifc,ie)
> end do
> end do
> facint2=sum
> return
> end
> c-----------------------------------------------------------------------
> subroutine out_csrmats(acsr,ia,ja,n,name9)
> real acsr(1)
> integer ia(1),ja(1)
> c
> character*9 name9
> character*9 s(16)
> c
> nnz = ia(n+1)-ia(1)
> c
> write(6,1) name9,n,nnz
> 1 format(/,'CSR Mat:',a9,3x,'n =',i5,3x,'nnz =',i5,/)
> c
> n16 = min(n,16)
> n29 = min(n,29)
> do i=1,n29
> call blank(s,9*16)
> n1 = ia(i)
> n2 = ia(i+1)-1
> do jj=n1,n2
> j = ja (jj)
> a = acsr(jj)
> if (a.ne.0..and.j.le.n16) write(s(j),9) a
> enddo
> write(6,16) (s(k),k=1,n16)
> enddo
> 9 format(f9.4)
> 16 format(16a9)
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine local_grad3(ur,us,ut,u,N,e,D,Dt)
> c Output: ur,us,ut Input:u,N,e,D,Dt
> real ur(0:N,0:N,0:N),us(0:N,0:N,0:N),ut(0:N,0:N,0:N)
> real u (0:N,0:N,0:N,1)
> real D (0:N,0:N),Dt(0:N,0:N)
> integer e
> c
> m1 = N+1
> m2 = m1*m1
> c
> call mxm(D ,m1,u(0,0,0,e),m1,ur,m2)
> do k=0,N
> call mxm(u(0,0,k,e),m1,Dt,m1,us(0,0,k),m1)
> enddo
> call mxm(u(0,0,0,e),m2,Dt,m1,ut,m1)
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine local_grad2(ur,us,u,N,e,D,Dt)
> c Output: ur,us Input:u,N,e,D,Dt
> real ur(0:N,0:N),us(0:N,0:N)
> real u (0:N,0:N,1)
> real D (0:N,0:N),Dt(0:N,0:N)
> integer e
> c
> m1 = N+1
> c
> call mxm(D ,m1,u(0,0,e),m1,ur,m1)
> call mxm(u(0,0,e),m1,Dt,m1,us,m1)
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine gradm1(ux,uy,uz,u)
> c
> c Compute gradient of T -- mesh 1 to mesh 1 (vel. to vel.)
> c
> include 'SIZE'
> include 'DXYZ'
> include 'GEOM'
> include 'INPUT'
> include 'TSTEP'
> c
> parameter (lxyz=lx1*ly1*lz1)
> real ux(lxyz,1),uy(lxyz,1),uz(lxyz,1),u(lxyz,1)
>
> common /ctmp1/ ur(lxyz),us(lxyz),ut(lxyz)
>
> integer e
>
> nxyz = nx1*ny1*nz1
> ntot = nxyz*nelt
>
> N = nx1-1
> do e=1,nelt
> if (if3d) then
> call local_grad3(ur,us,ut,u,N,e,dxm1,dxtm1)
> do i=1,lxyz
> ux(i,e) = jacmi(i,e)*(ur(i)*rxm1(i,1,1,e)
> $ + us(i)*sxm1(i,1,1,e)
> $ + ut(i)*txm1(i,1,1,e) )
> uy(i,e) = jacmi(i,e)*(ur(i)*rym1(i,1,1,e)
> $ + us(i)*sym1(i,1,1,e)
> $ + ut(i)*tym1(i,1,1,e) )
> uz(i,e) = jacmi(i,e)*(ur(i)*rzm1(i,1,1,e)
> $ + us(i)*szm1(i,1,1,e)
> $ + ut(i)*tzm1(i,1,1,e) )
> enddo
> else
> if (ifaxis) call setaxdy (ifrzer(e))
> call local_grad2(ur,us,u,N,e,dxm1,dytm1)
> do i=1,lxyz
> ux(i,e) =jacmi(i,e)*(ur(i)*rxm1(i,1,1,e)
> $ + us(i)*sxm1(i,1,1,e) )
> uy(i,e) =jacmi(i,e)*(ur(i)*rym1(i,1,1,e)
> $ + us(i)*sym1(i,1,1,e) )
> enddo
> endif
> enddo
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine comp_vort3(vort,work1,work2,u,v,w)
> c
> include 'SIZE'
> include 'TOTAL'
> c
> parameter(lt=lx1*ly1*lz1*lelv)
> real vort(lt,3),work1(1),work2(1),u(1),v(1),w(1)
> c
> ntot = nx1*ny1*nz1*nelv
> if (if3d) then
> c work1=dw/dy ; work2=dv/dz
> call dudxyz(work1,w,rym1,sym1,tym1,jacm1,1,2)
> call dudxyz(work2,v,rzm1,szm1,tzm1,jacm1,1,3)
> call sub3(vort(1,1),work1,work2,ntot)
> c work1=du/dz ; work2=dw/dx
> call dudxyz(work1,u,rzm1,szm1,tzm1,jacm1,1,3)
> call dudxyz(work2,w,rxm1,sxm1,txm1,jacm1,1,1)
> call sub3(vort(1,2),work1,work2,ntot)
> c work1=dv/dx ; work2=du/dy
> call dudxyz(work1,v,rxm1,sxm1,txm1,jacm1,1,1)
> call dudxyz(work2,u,rym1,sym1,tym1,jacm1,1,2)
> call sub3(vort(1,3),work1,work2,ntot)
> else
> c work1=dv/dx ; work2=du/dy
> call dudxyz(work1,v,rxm1,sxm1,txm1,jacm1,1,1)
> call dudxyz(work2,u,rym1,sym1,tym1,jacm1,1,2)
> call sub3(vort,work1,work2,ntot)
> endif
> c
> c Avg at bndry
> c
> ifielt = ifield
> ifield = 1
> if (if3d) then
> do idim=1,ndim
> call col2 (vort(1,idim),bm1,ntot)
> call dssum (vort(1,idim),nx1,ny1,nz1)
> call col2 (vort(1,idim),binvm1,ntot)
> enddo
> else
> call col2 (vort,bm1,ntot)
> call dssum (vort,nx1,ny1,nz1)
> call col2 (vort,binvm1,ntot)
> endif
> ifield = ifielt
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine surface_int(sint,sarea,a,ie,iface1)
> C
> include 'SIZE'
> include 'GEOM'
> include 'PARALLEL'
> include 'TOPOL'
> real a(lx1,ly1,lz1,1)
> c
> integer icalld
> save icalld
> data icalld/0/
> logical ifpf
> save ifpf
> c
> if (icalld.eq.0) then
> icalld=icalld+1
> if (skpdat(1,2).eq.nx1) then
> c write(6,*) 'In surface_int, using pf version of skpdat.'
> ifpf = .true.
> else
> c write(6,*) 'In surface_int, using std version of skpdat.'
> ifpf = .false.
> endif
> endif
> C
> sarea = 0.
> sint = 0.
> C
> call dsset(nx1,ny1,nz1)
> iface = eface1(iface1)
> c
> c Check skpdat (because of difference in pf vs. commercial
> version...arrghh)
> c
> if (ifpf) then
> c pf version
> js1 = skpdat(1,iface)
> jf1 = skpdat(2,iface)
> jskip1 = skpdat(3,iface)
> js2 = skpdat(4,iface)
> jf2 = skpdat(5,iface)
> jskip2 = skpdat(6,iface)
> else
> c std version
> js1 = skpdat(iface,1)
> jf1 = skpdat(iface,2)
> jskip1 = skpdat(iface,3)
> js2 = skpdat(iface,4)
> jf2 = skpdat(iface,5)
> jskip2 = skpdat(iface,6)
> endif
> C
> I = 0
> do 100 j2=js2,jf2,jskip2
> do 100 j1=js1,jf1,jskip1
> I = I+1
> sarea = sarea+area(i,1,iface1,ie)
> sint = sint +area(i,1,iface1,ie)*a(j1,j2,1,ie)
> 100 continue
> C
> return
> end
> c-----------------------------------------------------------------------
> subroutine surface_flux(dq,qx,qy,qz,ie,iface,w)
> C
> include 'SIZE'
> include 'GEOM'
> include 'INPUT'
> include 'PARALLEL'
> include 'TOPOL'
> parameter (l=lx1*ly1*lz1)
> real w(lx1,ly1,lz1),qx(l,1),qy(l,1),qz(l,1)
> c
> integer icalld
> save icalld
> data icalld/0/
> logical ifpf
> save ifpf
> c
> call dsset(nx1,ny1,nz1)
> if (icalld.eq.0) then
> icalld=icalld+1
> if (skpdat(1,2).eq.nx1) then
> c write(6,*) 'In surface_flux, using pf version of skpdat.'
> ifpf = .true.
> else
> c write(6,*) 'In surface_flux, using std version of skpdat.'
> ifpf = .false.
> endif
> endif
> C
> ifacepf = eface1(iface)
> c
> c Check skpdat (because of difference in pf vs. commercial
> version...arrghh)
> c
> if (ifpf) then
> c pf version
> js1 = skpdat(1,ifacepf)
> jf1 = skpdat(2,ifacepf)
> jskip1 = skpdat(3,ifacepf)
> js2 = skpdat(4,ifacepf)
> jf2 = skpdat(5,ifacepf)
> jskip2 = skpdat(6,ifacepf)
> else
> c std version
> js1 = skpdat(ifacepf,1)
> jf1 = skpdat(ifacepf,2)
> jskip1 = skpdat(ifacepf,3)
> js2 = skpdat(ifacepf,4)
> jf2 = skpdat(ifacepf,5)
> jskip2 = skpdat(ifacepf,6)
> endif
> C
> call faccl3 (w,qx(1,ie),unx(1,1,iface,ie),iface)
> call faddcl3(w,qy(1,ie),uny(1,1,iface,ie),iface)
> if (if3d)
> $call faddcl3(w,qz(1,ie),unz(1,1,iface,ie),iface)
> c
> dq = 0
> i = 0
> do 100 j2=js2,jf2,jskip2
> do 100 j1=js1,jf1,jskip1
> i = i+1
> dq = dq +area(i,1,iface,ie)*w(j1,j2,1)
> 100 continue
> C
> return
> end
> c-----------------------------------------------------------------------
> subroutine gaujordf(a,m,n,indr,indc,ipiv,ierr,rmult)
> C
> C Gauss-Jordan matrix inversion with full pivoting
> c
> c Num. Rec. p. 30, 2nd Ed., Fortran
> c
> C
> C a is an m x n matrix
> C rmult is a work array of dimension m
> C
> c
> real a(m,n),rmult(m)
> integer indr(m),indc(n),ipiv(n)
>
> c call outmat(a,m,n,'ab4',n)
> c do i=1,m
> c write(6,1) (a(i,j),j=1,n)
> c enddo
> c 1 format('mat: ',1p6e12.4)
>
> ierr = 0
> eps = 1.e-9
> call izero(ipiv,m)
>
> do k=1,m
> amx=0.
> do i=1,m ! Pivot search
> if (ipiv(i).ne.1) then
> do j=1,m
> if (ipiv(j).eq.0) then
> if (abs(a(i,j)).ge.amx) then
> amx = abs(a(i,j))
> ir = i
> jc = j
> endif
> elseif (ipiv(j).gt.1) then
> ierr = -ipiv(j)
> return
> endif
> enddo
> endif
> enddo
> ipiv(jc) = ipiv(jc) + 1
> c
> c Swap rows
> if (ir.ne.jc) then
> do j=1,n
> tmp = a(ir,j)
> a(ir,j) = a(jc,j)
> a(jc,j) = tmp
> enddo
> endif
> indr(k)=ir
> indc(k)=jc
> c write(6 ,*) k,' Piv:',jc,a(jc,jc)
> c write(28,*) k,' Piv:',jc,a(jc,jc)
> if (abs(a(jc,jc)).lt.eps) then
> write(6 ,*) 'small Gauss Jordan Piv:',jc,a(jc,jc)
> write(28,*) 'small Gauss Jordan Piv:',jc,a(jc,jc)
> ierr = jc
> call exitt
> return
> endif
> piv = 1./a(jc,jc)
> a(jc,jc)=1.
> do j=1,n
> a(jc,j) = a(jc,j)*piv
> enddo
> c
> do j=1,n
> work = a(jc,j)
> a(jc,j) = a(1 ,j)
> a(1 ,j) = work
> enddo
> do i=2,m
> rmult(i) = a(i,jc)
> a(i,jc) = 0.
> enddo
> c
> do j=1,n
> do i=2,m
> a(i,j) = a(i,j) - rmult(i)*a(1,j)
> enddo
> enddo
> c
> do j=1,n
> work = a(jc,j)
> a(jc,j) = a(1 ,j)
> a(1 ,j) = work
> enddo
> c
> c do i=1,m
> c if (i.ne.jc) then
> c rmult = a(i,jc)
> c a(i,jc) = 0.
> c do j=1,n
> c a(i,j) = a(i,j) - rmult*a(jc,j)
> c enddo
> c endif
> c enddo
> c
> enddo
> c
> c Unscramble matrix
> do j=m,1,-1
> if (indr(j).ne.indc(j)) then
> do i=1,m
> tmp=a(i,indr(j))
> a(i,indr(j))=a(i,indc(j))
> a(i,indc(j))=tmp
> enddo
> endif
> enddo
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine legendre_poly(L,x,N)
> c
> c Evaluate Legendre polynomials of degrees 0-N at point x
> c
> real L(0:N)
> c
> L(0) = 1.
> L(1) = x
> c
> do j=2,N
> L(j) = ( (2*j-1) * x * L(j-1) - (j-1) * L(j-2) ) / j
> enddo
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine build_new_filter(intv,zpts,nx,kut,wght,nid)
> c
> c This routing builds a 1D filter with a transfer function that
> c looks like:
> c
> c
> c ^
> c d_k |
> c | |
> c 1 |__________ _v_
> c | -_ c | \ wght
> c | \ ___
> c | | ^
> c 0 |-------------|---|>
> c
> c 0 c N k-->
> c
> c Where c := N-kut is the point below which d_k = 1.
> c
> c
> c
> c Here, nx = number of points
> c
> real intv(nx,nx),zpts(nx)
> c
> parameter (lm=40)
> parameter (lm2=lm*lm)
> real phi(lm2),pht(lm2),diag(lm2),rmult(lm),Lj(lm)
> integer indr(lm),indc(lm),ipiv(lm)
> c
> if (nx.gt.lm) then
> write(6,*) 'ABORT in build_new_filter:',nx,lm
> call exitt
> endif
> c
> kj = 0
> n = nx-1
> do j=1,nx
> z = zpts(j)
> call legendre_poly(Lj,z,n)
> kj = kj+1
> pht(kj) = Lj(1)
> kj = kj+1
> pht(kj) = Lj(2)
> do k=3,nx
> kj = kj+1
> pht(kj) = Lj(k)-Lj(k-2)
> enddo
> enddo
> call transpose (phi,nx,pht,nx)
> call copy (pht,phi,nx*nx)
> call gaujordf (pht,nx,nx,indr,indc,ipiv,ierr,rmult)
> c
> c Set up transfer function
> c
> call ident (diag,nx)
> c
> k0 = nx-kut
> do k=k0+1,nx
> kk = k+nx*(k-1)
> amp = wght*(k-k0)*(k-k0)/(kut*kut) ! quadratic growth
> diag(kk) = 1.-amp
> enddo
> c
> call mxm (diag,nx,pht,nx,intv,nx) ! -1
> call mxm (phi ,nx,intv,nx,pht,nx) ! V D V
> call copy (intv,pht,nx*nx)
> c
> do k=1,nx*nx
> pht(k) = 1.-diag(k)
> enddo
> np1 = nx+1
> if (nid.eq.0) then
> write(6,6) 'filt amp',(pht (k),k=1,nx*nx,np1)
> write(6,6) 'filt trn',(diag(k),k=1,nx*nx,np1)
> 6 format(a8,16f7.4,6(/,8x,16f7.4))
> endif
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine avg_all
> c
> c This routine computes running averages E(X),E(X^2),E(X*Y)
> c and outputs to avg*.fld*, rms*.fld*, and rm2*.fld* for all
> c fields.
> c
> c E denotes the expected value operator and X,Y two
> c real valued random variables.
> c
> c variances and covariances can be computed in a post-processing
> step:
> c
> c var(X) := E(X^X) - E(X)*E(X) c cov(X,Y) := E(X*Y) -
> E(X)*E(Y) c
> c Note: The E-operator is linear, in the sense that the expected
> c value is given by E(X) = 1/N * sum[ E(X)_i ], where E(X)_i
> c is the expected value of the sub-ensemble i (i=1...N).
> c
> include 'SIZE'
> include 'TOTAL'
> include 'AVG'
>
> logical ifverbose
> integer icalld
> save icalld
> data icalld /0/
>
> if (ax1.ne.lx1 .or. ay1.ne.ly1 .or. az1.ne.lz1) then
> if(nid.eq.0) write(6,*)
> $ 'ABORT: wrong size of ax1,ay1,az1 in avg_all(), check SIZE!'
> call exitt
> endif
> if (ax2.ne.lx2 .or. ay2.ne.ay2 .or. az2.ne.lz2) then
> if(nid.eq.0) write(6,*)
> $ 'ABORT: wrong size of ax2,ay2,az2 in avg_all(), check SIZE!'
> call exitt
> endif
>
> ntot = nx1*ny1*nz1*nelv
> ntott = nx1*ny1*nz1*nelt
> nto2 = nx2*ny2*nz2*nelv
>
> ! initialization
> if (icalld.eq.0) then
> icalld = icalld + 1
> atime = 0.
> timel = time
>
> call rzero(uavg,ntot)
> call rzero(vavg,ntot)
> call rzero(wavg,ntot)
> call rzero(pavg,nto2)
> do i = 1,ldimt
> call rzero(tavg(1,1,1,1,i),ntott)
> enddo
>
> call rzero(urms,ntot)
> call rzero(vrms,ntot)
> call rzero(wrms,ntot)
> call rzero(prms,nto2)
> do i = 1,ldimt
> call rzero(trms(1,1,1,1,i),ntott)
> enddo
>
> call rzero(vwms,ntot)
> call rzero(wums,ntot)
> call rzero(uvms,ntot)
> endif
>
> dtime = time - timel
> atime = atime + dtime
>
> ! dump freq
> iastep = param(68)
> if (iastep.eq.0) iastep=param(15) ! same as iostep
> if (iastep.eq.0) iastep=500
>
> ifverbose=.false.
> if (istep.le.10) ifverbose=.true.
> if (mod(istep,iastep).eq.0) ifverbose=.true.
>
> if (atime.ne.0..and.dtime.ne.0.) then
> if(nid.eq.0) write(6,*) 'Compute statistics ...'
> beta = dtime/atime
> alpha = 1.-beta
> ! compute averages E(X)
> call avg1 (uavg,vx,alpha,beta,ntot ,'um ',ifverbose)
> call avg1 (vavg,vy,alpha,beta,ntot ,'vm ',ifverbose)
> call avg1 (wavg,vz,alpha,beta,ntot ,'wm ',ifverbose)
> call avg1 (pavg,pr,alpha,beta,nto2 ,'prm ',ifverbose)
> call avg1 (tavg,t ,alpha,beta,ntott,'tm ',ifverbose)
> do i = 2,ldimt
> call avg1 (tavg(1,1,1,1,i),t(1,1,1,1,i),alpha,beta,
> & ntott,'psav',ifverbose)
> enddo
>
> ! compute averages E(X^2)
> call avg2 (urms,vx,alpha,beta,ntot ,'ums ',ifverbose)
> call avg2 (vrms,vy,alpha,beta,ntot ,'vms ',ifverbose)
> call avg2 (wrms,vz,alpha,beta,ntot ,'wms ',ifverbose)
> call avg2 (prms,pr,alpha,beta,nto2 ,'prms',ifverbose)
> call avg2 (trms,t ,alpha,beta,ntott,'tms ',ifverbose)
> do i = 2,ldimt
> call avg2 (trms(1,1,1,1,i),t(1,1,1,1,i),alpha,beta,
> & ntott,'psms',ifverbose)
> enddo
>
> ! compute averages E(X*Y) (for now just for the velocities)
> call avg3 (uvms,vx,vy,alpha,beta,ntot,'uvm ',ifverbose)
> call avg3 (vwms,vy,vz,alpha,beta,ntot,'vwm ',ifverbose)
> call avg3 (wums,vz,vx,alpha,beta,ntot,'wum ',ifverbose)
> endif
> c
> c-----------------------------------------------------------------------
> if ( (mod(istep,iastep).eq.0.and.istep.gt.1) .or.lastep.eq.1)
> then
>
> time_temp = time
> time = atime ! Output the duration of this avg
>
> call outpost2(uavg,vavg,wavg,pavg,tavg,ldimt,'avg')
> call outpost2(urms,vrms,wrms,prms,trms,ldimt,'rms')
> call outpost (uvms,vwms,wums,prms,trms, 'rm2')
>
> atime = 0.
> time = time_temp ! Restore clock
>
> endif
> c
> timel = time
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine avg1(avg,f,alpha,beta,n,name,ifverbose)
> include 'SIZE'
> include 'TSTEP'
> c
> real avg(n),f(n)
> character*4 name
> logical ifverbose
> c
> do k=1,n
> avg(k) = alpha*avg(k) + beta*f(k)
> enddo
> c
> if (ifverbose) then
> avgmax = glmax(avg,n)
> avgmin = glmin(avg,n)
> if (nid.eq.0) write(6,1) istep,time,avgmin,avgmax
> $ ,alpha,beta,name
> 1 format(i9,1p5e13.5,1x,a4,' av1mnx')
> endif
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine avg2(avg,f,alpha,beta,n,name,ifverbose)
> include 'SIZE'
> include 'TSTEP'
> c
> real avg(n),f(n)
> character*4 name
> logical ifverbose
> c
> do k=1,n
> avg(k) = alpha*avg(k) + beta*f(k)*f(k)
> enddo
> c
> if (ifverbose) then
> avgmax = glmax(avg,n)
> avgmin = glmin(avg,n)
> if (nid.eq.0) write(6,1) istep,time,avgmin,avgmax
> $ ,alpha,beta,name
> 1 format(i9,1p5e13.5,1x,a4,' av2mnx')
> endif
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine avg3(avg,f,g,alpha,beta,n,name,ifverbose)
> include 'SIZE'
> include 'TSTEP'
> c
> real avg(n),f(n),g(n)
> character*4 name
> logical ifverbose
> c
> do k=1,n
> avg(k) = alpha*avg(k) + beta*f(k)*g(k)
> enddo
> c
> if (ifverbose) then
> avgmax = glmax(avg,n)
> avgmin = glmin(avg,n)
> if (nid.eq.0) write(6,1) istep,time,avgmin,avgmax
> $ ,alpha,beta,name
> 1 format(i9,1p5e13.5,1x,a4,' av3mnx')
> endif
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine build_legend_transform(Lj,Ljt,zpts,nx)
> c
> real Lj(nx*nx),Ljt(nx*nx),zpts(nx)
> c
> parameter (lm=90)
> integer indr(lm),indc(lm),ipiv(lm)
> c
> if (nx.gt.lm) then
> write(6,*) 'ABORT in build_legend_transform:',nx,lm
> call exitt
> endif
> c
> j = 1
> n = nx-1
> do i=1,nx
> z = zpts(i)
> call legendre_poly(Lj(j),z,n) ! Return Lk(z), k=0,...,n
> j = j+nx
> enddo
> call transpose1(Lj,nx)
> c call outmat(Lj,nx,nx,'Lj ',n)
> c call exitt
> call gaujordf (Lj,nx,nx,indr,indc,ipiv,ierr,rmult)
> call transpose (Ljt,nx,Lj,nx)
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine local_err_est(err,u,nx,Lj,Ljt,uh,w,if3d)
> c
> c Local error estimates for u_e
> c
> include 'SIZE'
> real err(5,2),u(1),uh(nx,nx,nx),w(1),Lj(1),Ljt(1)
> logical if3d
> c
> call rzero(err,10)
> c
> nxyz = nx**ndim
> utot = vlsc2(u,u,nxyz)
> if (utot.eq.0) return
> c
> call tensr3(uh,nx,u,nx,Lj,Ljt,Ljt,w) ! Go to Legendre space
> c
> c
> c Get energy in low modes
> c
> m = nx-2
> c
> if (if3d) then
> amp2_l = 0.
> do k=1,m
> do j=1,m
> do i=1,m
> amp2_l = amp2_l + uh(i,j,k)**2
> enddo
> enddo
> enddo
> c
> c Energy in each spatial direction
> c
> amp2_t = 0
> do k=m+1,nx
> do j=1,m
> do i=1,m
> amp2_t = amp2_t + uh(i,j,k)**2
> enddo
> enddo
> enddo
> c
> amp2_s = 0
> do k=1,m
> do j=m+1,nx
> do i=1,m
> amp2_s = amp2_s + uh(i,j,k)**2
> enddo
> enddo
> enddo
> c
> amp2_r = 0
> do k=1,m
> do j=1,m
> do i=m+1,nx
> amp2_r = amp2_r + uh(i,j,k)**2
> enddo
> enddo
> enddo
> c
> amp2_h = 0
> do k=m+1,nx
> do j=m+1,nx
> do i=m+1,nx
> amp2_h = amp2_h + uh(i,j,k)**2
> enddo
> enddo
> enddo
> c
> etot = amp2_l + amp2_r + amp2_s + amp2_t + amp2_h
> c
> relr = amp2_r / (amp2_r + amp2_l)
> rels = amp2_s / (amp2_s + amp2_l)
> relt = amp2_t / (amp2_t + amp2_l)
> relh = (amp2_r + amp2_s + amp2_t + amp2_h) / etot
> c
> else
> k = 1
> amp2_l = 0.
> do j=1,m
> do i=1,m
> amp2_l = amp2_l + uh(i,j,k)**2
> enddo
> enddo
> if (amp2_l.eq.0) return
> c
> c Energy in each spatial direction
> c
> amp2_t = 0
> c
> amp2_s = 0
> do j=m+1,nx
> do i=1,m
> amp2_s = amp2_s + uh(i,j,k)**2
> enddo
> enddo
> c
> amp2_r = 0
> do j=1,m
> do i=m+1,nx
> amp2_r = amp2_r + uh(i,j,k)**2
> enddo
> enddo
> c
> amp2_h = 0
> do j=m+1,nx
> do i=m+1,nx
> amp2_h = amp2_h + uh(i,j,k)**2
> enddo
> enddo
> c
> etot = amp2_l + amp2_r + amp2_s + amp2_h
> c
> relr = amp2_r / (amp2_r + amp2_l)
> rels = amp2_s / (amp2_s + amp2_l)
> relt = 0
> relh = (amp2_r + amp2_s + amp2_h) / etot
> c
> endif
> c
> err (1,1) = sqrt(relr)
> err (2,1) = sqrt(rels)
> if (if3d) err (3,1) = sqrt(relt)
> err (4,1) = sqrt(relh)
> err (5,1) = sqrt(etot)
> c
> err (1,2) = sqrt(amp2_r)
> err (2,2) = sqrt(amp2_s)
> if (if3d) err (3,2) = sqrt(amp2_t)
> err (4,2) = sqrt(amp2_h)
> err (5,2) = sqrt(utot)
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine transpose1(a,n)
> real a(n,n)
> c
> do j=1,n
> do i=j+1,n
> ta = a(i,j)
> a(i,j) = a(j,i)
> a(j,i) = ta
> enddo
> enddo
> return
> end
> c-----------------------------------------------------------------------
> subroutine get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
> integer ex,ey,ez,eg
> c
> nelxy = nelx*nely
> c
> ez = 1 + (eg-1)/nelxy
> ey = mod1 (eg,nelxy)
> ey = 1 + (ey-1)/nelx
> ex = mod1 (eg,nelx)
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine dump_header2d(excode,nx,ny,nlx,nly,ierr)
>
> include 'SIZE'
> include 'TOTAL'
>
> character*2 excode(15)
>
> real*4 test_pattern
>
> character*1 fhdfle1(132)
> character*132 fhdfle
> equivalence (fhdfle,fhdfle1)
>
> jstep = istep
> if (jstep.gt.10000) jstep = jstep / 10
> if (jstep.gt.10000) jstep = jstep / 10
> if (jstep.gt.10000) jstep = jstep / 10
> if (jstep.gt.10000) jstep = jstep / 10
> if (jstep.gt.10000) jstep = jstep / 10
>
> nlxy = nlx*nly
> nzz = 1
>
> c write(6,'(4i4,1PE14.7,i5,1x,15a2,1x,a12)')
> c $ nlxy,nx,ny,nzz,TIME,jstep,(excode(i),i=1,15),
> c $ 'NELT,NX,NY,N'
> c
> p66 = 0.
> ierr= 0
> IF (p66.EQ.1.0) THEN
> C unformatted i/o
> WRITE(24)
> $ nlxy,nx,ny,nzz,TIME,jstep,(excode(i),i=1,15)
> ELSEIF (p66.EQ.3.0) THEN
> C formatted i/o to header file
> WRITE(27,'(4I4,1pe14.7,I5,1X,15A2,1X,A12)')
> $ nlxy,nx,ny,nzz,TIME,jstep,(excode(i),i=1,15),
> $ 'NELT,NX,NY,N'
> ELSEIF (p66.eq.4.0) THEN
> C formatted i/o to header file
> WRITE(fhdfle,'(4I4,1pe14.7,I5,1X,15A2,1X,A12)')
> $ nlxy,nx,ny,nzz,TIME,jstep,(excode(i),i=1,15),
> $ ' 4 NELT,NX,NY,N'
> call byte_write(fhdfle,20,ierr)
> ELSEIF (p66.eq.5.0) THEN
> C formatted i/o to header file
> WRITE(fhdfle,'(4I4,1pe14.7,I5,1X,15A2,1X,A12)')
> $ nlxy,nx,ny,nzz,TIME,jstep,(excode(i),i=1,15),
> $ ' 8 NELT,NX,NY,N'
> call byte_write(fhdfle,20,ierr)
> ELSE
> C formatted i/o
> WRITE(24,'(4I4,1pe14.7,I5,1X,15A2,1X,A12)')
> $ nlxy,nx,ny,nzz,TIME,jstep,(excode(i),i=1,15),
> $ 'NELT,NX,NY,N'
> ENDIF
> C cdrror is a dummy cerror value for now.
> CDRROR=0.0
> IF (p66.EQ.1.0) THEN
> C unformatted i/o
> WRITE(24)(CDRROR,I=1,nlxy)
> ELSEIF (p66.eq.3. .or. p66.eq.4.0) then
> C write byte-ordering test pattern to byte file...
> test_pattern = 6.54321
> call byte_write(test_pattern,1,ierr)
> ELSEIF (p66.eq.5.) then
> test_pattern8 = 6.54321
> call byte_write(test_pattern8,2,ierr)
> ELSE
> C formatted i/o
> WRITE(24,'(6G11.4)')(CDRROR,I=1,nlxy)
> ENDIF
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine outfld2d_p(u,v,w,nx,ny,nlx,nly,name,ifld,jid,npido,ir)
>
> include 'SIZE'
> include 'TOTAL'
>
> integer icalld
> save icalld
> data icalld /0/
>
> real u(nx*ny*nlx*nly)
> real v(nx*ny*nlx*nly)
> real w(nx*ny*nlx*nly)
> character*4 name
>
> character*2 excode(15)
> character*12 fm
> character*20 outfile
>
> character*1 slash,dot
> save slash,dot
> data slash,dot / '/' , '.' /
>
> icalld = icalld+1
>
> call blank(excode,30)
> excode(4) = 'U '
> excode(5) = ' '
> excode(6) = 'T '
> nthings = 3
> ir = 0 !error code for dump_header2d
>
> call blank(outfile,20)
> if (npido.lt.100) then
> if (ifld.lt.100) then
> write(outfile,22) jid,slash,name,ifld
> 22 format('B',i2.2,a1,a4,'.fld',i2.2)
> elseif (ifld.lt.1000) then
> write(outfile,23) jid,slash,name,ifld
> 23 format('B',i2.2,a1,a4,'.fld',i3)
> elseif (ifld.lt.10000) then
> write(outfile,24) jid,slash,name,ifld
> 24 format('B',i2.2,a1,a4,'.fld',i4)
> elseif (ifld.lt.100000) then
> write(outfile,25) jid,slash,name,ifld
> 25 format('B',i2.2,a1,a4,'.fld',i5)
> elseif (ifld.lt.1000000) then
> write(outfile,26) jid,slash,name,ifld
> 26 format('B',i2.2,a1,a4,'.fld',i6)
> endif
> else
> if (ifld.lt.100) then
> write(outfile,32) jid,slash,name,ifld
> 32 format('B',i3.3,a1,a4,'.fld',i2.2)
> elseif (ifld.lt.1000) then
> write(outfile,33) jid,slash,name,ifld
> 33 format('B',i3.3,a1,a4,'.fld',i3)
> elseif (ifld.lt.10000) then
> write(outfile,34) jid,slash,name,ifld
> 34 format('B',i3.3,a1,a4,'.fld',i4)
> elseif (ifld.lt.100000) then
> write(outfile,35) jid,slash,name,ifld
> 35 format('B',i3.3,a1,a4,'.fld',i5)
> elseif (ifld.lt.1000000) then
> write(outfile,36) jid,slash,name,ifld
> 36 format('B',i3.3,a1,a4,'.fld',i6)
> endif
> endif
>
> if (icalld.le.4) write(6,*) nid,outfile,' OPEN',nlx,nly
> open(unit=24,file=outfile,status='unknown')
> call dump_header2d(excode,nx,ny,nlx,nly,ir)
>
> n = nx*ny*nlx*nly
> write(fm,10) nthings
> write(24,fm) (u(i),v(i),w(i),i=1,n)
> 10 format('(1p',i1,'e14.6)')
> close(24)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine outfld2d(u,v,w,nx,ny,nlx,nly,name,ifld)
>
> include 'SIZE'
> include 'TOTAL'
>
> real u(nx*ny*nlx*nly)
> real v(nx*ny*nlx*nly)
> real w(nx*ny*nlx*nly)
> character*3 name
>
> character*2 excode(15)
> character*12 fm
> character*20 outfile
>
> c if (istep.le.10) write(6,*) nid,' in out2d:',iz
>
> call blank(excode,30)
> c
> c excode(1) = 'X '
> c excode(2) = 'Y '
> c excode(3) = 'U '
> c excode(4) = 'V '
> c excode(5) = 'P '
> c excode(6) = 'T '
> c
> excode(4) = 'U '
> excode(5) = ' '
> excode(6) = 'T '
> nthings = 3
>
> ierr = 0
> if (nid.eq.0) then
> call blank(outfile,20)
> if (ifld.lt.100) then
> write(outfile,2) name,ifld
> 2 format(a3,'2d.fld',i2.2)
> elseif (ifld.lt.1000) then
> write(outfile,3) name,ifld
> 3 format(a3,'2d.fld',i3)
> elseif (ifld.lt.10000) then
> write(outfile,4) name,ifld
> 4 format(a3,'2d.fld',i4)
> elseif (ifld.lt.100000) then
> write(outfile,5) name,ifld
> 5 format(a3,'2d.fld',i5)
> elseif (ifld.lt.1000000) then
> write(outfile,6) name,ifld
> 6 format(a3,'2d.fld',i6)
> endif
> open(unit=24,file=outfile,status='unknown')
> call dump_header2d(excode,nx,ny,nlx,nly,ierr)
>
> n = nx*ny*nlx*nly
> write(fm,10) nthings
> c write(6,*) fm
> c call exitt
> write(24,fm) (u(i),v(i),w(i),i=1,n)
> 10 format('(1p',i1,'e14.6)')
> c 10 format('''(1p',i1,'e15.7)''')
> c 10 format(1p7e15.7)
> c
> close(24)
> endif
> call err_chk(ierr,'Error using byte_write,outfld2d. Abort.$')
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine planar_average_z(ua,u,w1,w2)
> c
> c Compute r-s planar average of quantity u()
> c
> include 'SIZE'
> include 'GEOM'
> include 'PARALLEL'
> include 'WZ'
> include 'ZPER'
> c
> real ua(nz1,nelz),u(nx1*ny1,nz1,nelv),w1(nz1,nelz),w2(nz1,nelz)
> integer e,eg,ez
> c
> melxy = nelx*nely
> c
> nz = nz1*nelz
> call rzero(ua,nz)
> call rzero(w1,nz)
> c
> do e=1,nelt
> c
> eg = lglel(e)
> ez = 1 + (eg-1)/melxy
> c
> do k=1,nz1
> do i=1,nx1*ny1
> zz = (1.-zgm1(k,3))/2. ! = 1 for k=1, = 0 for k=nz1
> aa = zz*area(i,1,5,e) + (1-zz)*area(i,1,6,e) ! wgtd
> jacobian
> w1(k,ez) = w1(k,ez) + aa
> ua(k,ez) = ua(k,ez) + aa*u(i,k,e)
> enddo
> enddo
> enddo
> c
> call gop(ua,w2,'+ ',nz)
> call gop(w1,w2,'+ ',nz)
> c
> do i=1,nz
> ua(i,1) = ua(i,1) / w1(i,1) ! Normalize
> enddo
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine drgtrq(dgtq,xm0,ym0,zm0,sij,pm1,visc,f,e)
> c
> INCLUDE 'SIZE'
> INCLUDE 'GEOM'
> INCLUDE 'INPUT'
> INCLUDE 'TOPOL'
> INCLUDE 'TSTEP'
> c
> real dgtq(3,4)
> real xm0 (lx1,ly1,lz1,lelt)
> real ym0 (lx1,ly1,lz1,lelt)
> real zm0 (lx1,ly1,lz1,lelt)
> real sij (lx1,ly1,lz1,3*ldim-3,lelv)
> real pm1 (lx1,ly1,lz1,lelv)
> real visc(lx1,ly1,lz1,lelv)
> c
> real dg(3,2)
> c
> integer f,e,pf
> real n1,n2,n3
> c
> call dsset(nx1,ny1,nz1) ! set up counters
> pf = eface1(f) ! convert from preproc. notation
> js1 = skpdat(1,pf)
> jf1 = skpdat(2,pf)
> jskip1 = skpdat(3,pf)
> js2 = skpdat(4,pf)
> jf2 = skpdat(5,pf)
> jskip2 = skpdat(6,pf)
> C
> call rzero(dgtq,12)
> c
> if (if3d.or.ifaxis) then
> i = 0
> a = 0
> do j2=js2,jf2,jskip2
> do j1=js1,jf1,jskip1
> i = i+1
> n1 = unx(i,1,f,e)*area(i,1,f,e)
> n2 = uny(i,1,f,e)*area(i,1,f,e)
> n3 = unz(i,1,f,e)*area(i,1,f,e)
> a = a + area(i,1,f,e)
> c
> v = visc(j1,j2,1,e)
> c
> s11 = sij(j1,j2,1,1,e)
> s21 = sij(j1,j2,1,4,e)
> s31 = sij(j1,j2,1,6,e)
> c
> s12 = sij(j1,j2,1,4,e)
> s22 = sij(j1,j2,1,2,e)
> s32 = sij(j1,j2,1,5,e)
> c
> s13 = sij(j1,j2,1,6,e)
> s23 = sij(j1,j2,1,5,e)
> s33 = sij(j1,j2,1,3,e)
> c
> dg(1,1) = pm1(j1,j2,1,e)*n1 ! pressure drag
> dg(2,1) = pm1(j1,j2,1,e)*n2
> dg(3,1) = pm1(j1,j2,1,e)*n3
> c
> dg(1,2) = -v*(s11*n1 + s12*n2 + s13*n3) ! viscous drag
> dg(2,2) = -v*(s21*n1 + s22*n2 + s23*n3)
> dg(3,2) = -v*(s31*n1 + s32*n2 + s33*n3)
> c
> r1 = xm0(j1,j2,1,e)
> r2 = ym0(j1,j2,1,e)
> r3 = zm0(j1,j2,1,e)
> c
> do l=1,2
> do k=1,3
> dgtq(k,l) = dgtq(k,l) + dg(k,l)
> enddo
> enddo
> c
> dgtq(1,3) = dgtq(1,3) + (r2*dg(3,1)-r3*dg(2,1)) ! pressure
> dgtq(2,3) = dgtq(2,3) + (r3*dg(1,1)-r1*dg(3,1)) ! torque
> dgtq(3,3) = dgtq(3,3) + (r1*dg(2,1)-r2*dg(1,1))
> c
> dgtq(1,4) = dgtq(1,4) + (r2*dg(3,2)-r3*dg(2,2)) ! viscous
> dgtq(2,4) = dgtq(2,4) + (r3*dg(1,2)-r1*dg(3,2)) ! torque
> dgtq(3,4) = dgtq(3,4) + (r1*dg(2,2)-r2*dg(1,2))
> enddo
> enddo
>
> else ! 2D
>
> i = 0
> a = 0
> do j2=js2,jf2,jskip2
> do j1=js1,jf1,jskip1
> i = i+1
> n1 = unx(i,1,f,e)*area(i,1,f,e)
> n2 = uny(i,1,f,e)*area(i,1,f,e)
> a = a + area(i,1,f,e)
> v = visc(j1,j2,1,e)
>
> s11 = sij(j1,j2,1,1,e)
> s12 = sij(j1,j2,1,3,e)
> s21 = sij(j1,j2,1,3,e)
> s22 = sij(j1,j2,1,2,e)
>
> dg(1,1) = pm1(j1,j2,1,e)*n1 ! pressure drag
> dg(2,1) = pm1(j1,j2,1,e)*n2
> dg(3,1) = 0
>
> dg(1,2) = -v*(s11*n1 + s12*n2) ! viscous drag
> dg(2,2) = -v*(s21*n1 + s22*n2)
> dg(3,2) = 0.
>
> r1 = xm0(j1,j2,1,e)
> r2 = ym0(j1,j2,1,e)
> r3 = 0.
>
> do l=1,2
> do k=1,3
> dgtq(k,l) = dgtq(k,l) + dg(k,l)
> enddo
> enddo
>
> dgtq(1,3) = 0! dgtq(1,3) + (r2*dg(3,1)-r3*dg(2,1)) ! pressure
> dgtq(2,3) = 0! dgtq(2,3) + (r3*dg(1,1)-r1*dg(3,1)) ! torque
> dgtq(3,3) = dgtq(3,3) + (r1*dg(2,1)-r2*dg(1,1))
>
> dgtq(1,4) = 0! dgtq(1,4) + (r2*dg(3,2)-r3*dg(2,2)) ! viscous
> dgtq(2,4) = 0! dgtq(2,4) + (r3*dg(1,2)-r1*dg(3,2)) ! torque
> dgtq(3,4) = dgtq(3,4) + (r1*dg(2,2)-r2*dg(1,2))
> enddo
> enddo
> endif
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine torque_calc(scale,x0,ifdout,iftout)
> c
> c Compute torque about point x0
> c
> c Scale is a user-supplied multiplier so that results may be
> c scaled to any convenient non-dimensionalization.
> c
> c
> INCLUDE 'SIZE'
> INCLUDE 'TOTAL'
>
> common /cvflow_r/ flow_rate,base_flow,domain_length,xsec
> $ , scale_vf(3)
>
> c
> real x0(3),w1(0:maxobj)
> logical ifdout,iftout
> c
> common /scrns/ sij (lx1*ly1*lz1*6*lelv)
> common /scrcg/ pm1 (lx1,ly1,lz1,lelv)
> common /scrsf/ xm0(lx1,ly1,lz1,lelt)
> $, ym0(lx1,ly1,lz1,lelt)
> $, zm0(lx1,ly1,lz1,lelt)
> c
> parameter (lr=lx1*ly1*lz1)
> common /scruz/ ur(lr),us(lr),ut(lr)
> $ , vr(lr),vs(lr),vt(lr)
> $ , wr(lr),ws(lr),wt(lr)
> c
> common /ctorq/ dragx(0:maxobj),dragpx(0:maxobj),dragvx(0:maxobj)
> $ , dragy(0:maxobj),dragpy(0:maxobj),dragvy(0:maxobj)
> $ , dragz(0:maxobj),dragpz(0:maxobj),dragvz(0:maxobj)
> c
> $ , torqx(0:maxobj),torqpx(0:maxobj),torqvx(0:maxobj)
> $ , torqy(0:maxobj),torqpy(0:maxobj),torqvy(0:maxobj)
> $ , torqz(0:maxobj),torqpz(0:maxobj),torqvz(0:maxobj)
> c
> $ , dpdx_mean,dpdy_mean,dpdz_mean
> $ , dgtq(3,4)
> c
> c
> n = nx1*ny1*nz1*nelv
> c
> call mappr(pm1,pr,xm0,ym0) ! map pressure onto Mesh 1
> c
> c Add mean_pressure_gradient.X to p:
>
> if (param(55).ne.0) then
> dpdx_mean = -scale_vf(1)
> dpdy_mean = -scale_vf(2)
> dpdz_mean = -scale_vf(3)
> endif
>
> call add2s2(pm1,xm1,dpdx_mean,n) ! Doesn't work if object is cut
> by
> call add2s2(pm1,ym1,dpdy_mean,n) ! periodicboundary. In this
> case,
> call add2s2(pm1,zm1,dpdz_mean,n) ! set ._mean=0 and compensate
> in
> c
> c Compute sij
> c
> nij = 3
> if (if3d.or.ifaxis) nij=6
> call comp_sij(sij,nij,vx,vy,vz,ur,us,ut,vr,vs,vt,wr,ws,wt)
> c
> c
> c Fill up viscous array w/ default
> c
> if (istep.lt.1) call cfill(vdiff,param(2),n)
> c
> call cadd2(xm0,xm1,-x0(1),n)
> call cadd2(ym0,ym1,-x0(2),n)
> call cadd2(zm0,zm1,-x0(3),n)
> c
> x1min=glmin(xm0(1,1,1,1),n)
> x2min=glmin(ym0(1,1,1,1),n)
> x3min=glmin(zm0(1,1,1,1),n)
> c
> x1max=glmax(xm0(1,1,1,1),n)
> x2max=glmax(ym0(1,1,1,1),n)
> x3max=glmax(zm0(1,1,1,1),n)
> c
> do i=0,maxobj
> dragpx(i) = 0 ! BIG CODE :}
> dragvx(i) = 0
> dragx (i) = 0
> dragpy(i) = 0
> dragvy(i) = 0
> dragy (i) = 0
> dragpz(i) = 0
> dragvz(i) = 0
> dragz (i) = 0
> torqpx(i) = 0
> torqvx(i) = 0
> torqx (i) = 0
> torqpy(i) = 0
> torqvy(i) = 0
> torqy (i) = 0
> torqpz(i) = 0
> torqvz(i) = 0
> torqz (i) = 0
> enddo
> c
> c
> nobj = 0
> do ii=1,nhis
> if (hcode(10,ii).EQ.'I') then
> iobj = lochis(1,ii)
> memtot = nmember(iobj)
> nobj = max(iobj,nobj)
> c
> if (hcode(1,ii).ne.' ' .or. hcode(2,ii).ne.' ' .or.
> $ hcode(3,ii).ne.' ' ) then
> ifield = 1
> c
> c Compute drag for this object
> c
> do mem=1,memtot
> ieg = object(iobj,mem,1)
> ifc = object(iobj,mem,2)
> if (gllnid(ieg).eq.nid) then ! this processor has a
> contribution
> ie = gllel(ieg)
> call drgtrq(dgtq,xm0,ym0,zm0,sij,pm1,vdiff,ifc,ie)
> c
> call cmult(dgtq,scale,12)
> c
> dragpx(iobj) = dragpx(iobj) + dgtq(1,1) ! pressure
> dragpy(iobj) = dragpy(iobj) + dgtq(2,1)
> dragpz(iobj) = dragpz(iobj) + dgtq(3,1)
> c
> dragvx(iobj) = dragvx(iobj) + dgtq(1,2) ! viscous
> dragvy(iobj) = dragvy(iobj) + dgtq(2,2)
> dragvz(iobj) = dragvz(iobj) + dgtq(3,2)
> c
> torqpx(iobj) = torqpx(iobj) + dgtq(1,3) ! pressure
> torqpy(iobj) = torqpy(iobj) + dgtq(2,3)
> torqpz(iobj) = torqpz(iobj) + dgtq(3,3)
> c
> torqvx(iobj) = torqvx(iobj) + dgtq(1,4) ! viscous
> torqvy(iobj) = torqvy(iobj) + dgtq(2,4)
> torqvz(iobj) = torqvz(iobj) + dgtq(3,4)
> c
> endif
> enddo
> endif
> endif
> enddo
> c
> c Sum contributions from all processors
> c
> call gop(dragpx,w1,'+ ',maxobj+1)
> call gop(dragpy,w1,'+ ',maxobj+1)
> call gop(dragpz,w1,'+ ',maxobj+1)
> call gop(dragvx,w1,'+ ',maxobj+1)
> call gop(dragvy,w1,'+ ',maxobj+1)
> call gop(dragvz,w1,'+ ',maxobj+1)
> c
> call gop(torqpx,w1,'+ ',maxobj+1)
> call gop(torqpy,w1,'+ ',maxobj+1)
> call gop(torqpz,w1,'+ ',maxobj+1)
> call gop(torqvx,w1,'+ ',maxobj+1)
> call gop(torqvy,w1,'+ ',maxobj+1)
> call gop(torqvz,w1,'+ ',maxobj+1)
> c
> nobj = iglmax(nobj,1)
> c
> do i=1,nobj
> dragx(i) = dragpx(i) + dragvx(i)
> dragy(i) = dragpy(i) + dragvy(i)
> dragz(i) = dragpz(i) + dragvz(i)
> c
> torqx(i) = torqpx(i) + torqvx(i)
> torqy(i) = torqpy(i) + torqvy(i)
> torqz(i) = torqpz(i) + torqvz(i)
> c
> dragpx(0) = dragpx (0) + dragpx (i)
> dragvx(0) = dragvx (0) + dragvx (i)
> dragx (0) = dragx (0) + dragx (i)
> c
> dragpy(0) = dragpy (0) + dragpy (i)
> dragvy(0) = dragvy (0) + dragvy (i)
> dragy (0) = dragy (0) + dragy (i)
> c
> dragpz(0) = dragpz (0) + dragpz (i)
> dragvz(0) = dragvz (0) + dragvz (i)
> dragz (0) = dragz (0) + dragz (i)
> c
> torqpx(0) = torqpx (0) + torqpx (i)
> torqvx(0) = torqvx (0) + torqvx (i)
> torqx (0) = torqx (0) + torqx (i)
> c
> torqpy(0) = torqpy (0) + torqpy (i)
> torqvy(0) = torqvy (0) + torqvy (i)
> torqy (0) = torqy (0) + torqy (i)
> c
> torqpz(0) = torqpz (0) + torqpz (i)
> torqvz(0) = torqvz (0) + torqvz (i)
> torqz (0) = torqz (0) + torqz (i)
> c
> enddo
> c
> i0 = 0
> if (nobj.le.1) i0 = 1 ! one output for single-object case
> c
> do i=i0,nobj
> if (nid.eq.0) then
> if (if3d.or.ifaxis) then
> if (ifdout) then
> write(6,6)
> istep,time,dragx(i),dragpx(i),dragvx(i),i,'dragx'
> write(6,6)
> istep,time,dragy(i),dragpy(i),dragvy(i),i,'dragy'
> write(6,6)
> istep,time,dragz(i),dragpz(i),dragvz(i),i,'dragz'
> endif
> if (iftout) then
> write(6,6)
> istep,time,torqx(i),torqpx(i),torqvx(i),i,'torqx'
> write(6,6)
> istep,time,torqy(i),torqpy(i),torqvy(i),i,'torqy'
> write(6,6)
> istep,time,torqz(i),torqpz(i),torqvz(i),i,'torqz'
> endif
> else
> if (ifdout) then
> write(6,6)
> istep,time,dragx(i),dragpx(i),dragvx(i),i,'dragx'
> write(6,6)
> istep,time,dragy(i),dragpy(i),dragvy(i),i,'dragy'
> endif
> if (iftout) then
> write(6,6)
> istep,time,torqz(i),torqpz(i),torqvz(i),i,'torqz'
> endif
> endif
> endif
> 6 format(i8,1p4e19.11,2x,i5,a5)
> enddo
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine comp_sij(sij,nij,u,v,w,ur,us,ut,vr,vs,vt,wr,ws,wt)
> c du_i du_j
> c Compute the stress tensor S_ij := ---- + ----
> c du_j du_i
> c
> include 'SIZE'
> include 'TOTAL'
> c
> integer e
> c
> real sij(lx1*ly1*lz1,nij,lelv)
> real u (lx1*ly1*lz1,lelv)
> real v (lx1*ly1*lz1,lelv)
> real w (lx1*ly1*lz1,lelv)
> real ur (1) , us (1) , ut (1)
> $ , vr (1) , vs (1) , vt (1)
> $ , wr (1) , ws (1) , wt (1)
>
> real j ! Inverse Jacobian
>
> n = nx1-1 ! Polynomial degree
> nxyz = nx1*ny1*nz1
>
> if (if3d) then ! 3D CASE
> do e=1,nelv
> call local_grad3(ur,us,ut,u,N,e,dxm1,dxtm1)
> call local_grad3(vr,vs,vt,v,N,e,dxm1,dxtm1)
> call local_grad3(wr,ws,wt,w,N,e,dxm1,dxtm1)
>
> do i=1,nxyz
>
> j = jacmi(i,e)
>
> sij(i,1,e) = j* ! du/dx + du/dx
> $
> 2*(ur(i)*rxm1(i,1,1,e)+us(i)*sxm1(i,1,1,e)+ut(i)*txm1(i,1,1,e))
>
> sij(i,2,e) = j* ! dv/dy + dv/dy
> $
> 2*(vr(i)*rym1(i,1,1,e)+vs(i)*sym1(i,1,1,e)+vt(i)*tym1(i,1,1,e))
>
> sij(i,3,e) = j* ! dw/dz + dw/dz
> $
> 2*(wr(i)*rzm1(i,1,1,e)+ws(i)*szm1(i,1,1,e)+wt(i)*tzm1(i,1,1,e))
>
> sij(i,4,e) = j* ! du/dy + dv/dx
> $ (ur(i)*rym1(i,1,1,e)+us(i)*sym1(i,1,1,e)+ut(i)*tym1(i,1,1,e) +
> $ vr(i)*rxm1(i,1,1,e)+vs(i)*sxm1(i,1,1,e)+vt(i)*txm1(i,1,1,e) )
>
> sij(i,5,e) = j* ! dv/dz + dw/dy
> $ (wr(i)*rym1(i,1,1,e)+ws(i)*sym1(i,1,1,e)+wt(i)*tym1(i,1,1,e) +
> $ vr(i)*rzm1(i,1,1,e)+vs(i)*szm1(i,1,1,e)+vt(i)*tzm1(i,1,1,e) )
>
> sij(i,6,e) = j* ! du/dz + dw/dx
> $ (ur(i)*rzm1(i,1,1,e)+us(i)*szm1(i,1,1,e)+ut(i)*tzm1(i,1,1,e) +
> $ wr(i)*rxm1(i,1,1,e)+ws(i)*sxm1(i,1,1,e)+wt(i)*txm1(i,1,1,e) )
>
> enddo
> enddo
>
> elseif (ifaxis) then ! AXISYMMETRIC CASE
>
> c
> c Notation: ( 2 x Acheson, p. 353)
> c Cylindrical
> c Nek5k Coordinates
> c
> c x z
> c y r
> c z theta
> c
>
> do e=1,nelv
> call setaxdy ( ifrzer(e) ) ! change dytm1 if on-axis
> call local_grad2(ur,us,u,N,e,dxm1,dytm1)
> call local_grad2(vr,vs,v,N,e,dxm1,dytm1)
> call local_grad2(wr,ws,w,N,e,dxm1,dytm1)
>
> do i=1,nxyz
> j = jacmi(i,e)
> r = ym1(i,1,1,e) ! Cyl.
> Coord:
>
> sij(i,1,e) = j* ! du/dx + du/dx ! e_zz
> $ 2*(ur(i)*rxm1(i,1,1,e)+us(i)*sxm1(i,1,1,e))
>
> sij(i,2,e) = j* ! dv/dy + dv/dy ! e_rr
> $ 2*(vr(i)*rym1(i,1,1,e)+vs(i)*sym1(i,1,1,e))
>
> if (r.gt.0) then ! e_@@
> sij(i,3,e) = v(i,e)/r ! v / r
> else
> sij(i,3,e) = j* ! L'Hopital's rule: e_@@ = dv/dr
> $ 2*(vr(i)*rym1(i,1,1,e)+vs(i)*sym1(i,1,1,e))
> endif
>
> sij(i,4,e) = j* ! du/dy + dv/dx ! e_zr
> $ ( ur(i)*rym1(i,1,1,e)+us(i)*sym1(i,1,1,e) +
> $ vr(i)*rxm1(i,1,1,e)+vs(i)*sxm1(i,1,1,e) )
>
> if (yyyr.gt.0) then ! e_r@
> sij(i,5,e) = j* ! dw/dy
> $ ( wr(i)*rym1(i,1,1,e)+ws(i)*sym1(i,1,1,e) )
> $ - w(i,e) / r
> else
> sij(i,5,e) = 0
> endif
>
> sij(i,6,e) = j* ! dw/dx ! e_ at z
> $ ( wr(i)*rxm1(i,1,1,e)+ws(i)*sxm1(i,1,1,e) )
>
> enddo
> enddo
>
> else ! 2D CASE
>
> do e=1,nelv
> call local_grad2(ur,us,u,N,e,dxm1,dxtm1)
> call local_grad2(vr,vs,v,N,e,dxm1,dxtm1)
>
> do i=1,nxyz
> j = jacmi(i,e)
>
> sij(i,1,e) = j* ! du/dx + du/dx
> $ 2*(ur(i)*rxm1(i,1,1,e)+us(i)*sxm1(i,1,1,e))
>
> sij(i,2,e) = j* ! dv/dy + dv/dy
> $ 2*(vr(i)*rym1(i,1,1,e)+vs(i)*sym1(i,1,1,e))
>
> sij(i,3,e) = j* ! du/dy + dv/dx
> $ (ur(i)*rym1(i,1,1,e)+us(i)*sym1(i,1,1,e) +
> $ vr(i)*rxm1(i,1,1,e)+vs(i)*sxm1(i,1,1,e) )
>
> enddo
> enddo
> endif
> return
> end
> c-----------------------------------------------------------------------
> subroutine y_slice (ua,u,w1,w2)
> c
> c Extract a y slice of quantity u() - assumes global tens.prod.
> c
> include 'SIZE'
> include 'GEOM'
> include 'PARALLEL'
> include 'WZ'
> include 'ZPER'
> c
> real ua(nx1,nz1,nelx,nelz),u (nx1,ny1,nz1,nelv)
> $ ,w1(nx1,nz1,nelx,nelz),w2(nx1,nz1,nelx,nelz)
> integer e,eg,ex,ey,ez
> real dy2
> c
> mxz = nelx*nelz*nx1*nz1
> call rzero(ua,mxz)
> c
> do e=1,nelt
> c
> eg = lglel(e)
> call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
>
> j = 1
> if (ey.eq.1) then
> do k=1,nz1
> do i=1,nx1
> ua(i,k,ex,ez) = u(i,j,k,e)
> enddo
> enddo
> endif
> enddo
>
> call gop(ua,w2,'+ ',mxz)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine z_slice_g (uz,u,w1,kz,ezi,nx,ny,nz,nlxy)
> c
> c Extract a z slice of quantity u() c
> c ASSUMES data is in a global tensor-product structure, nlxy x
> nelz,
> c as would be produced by n2to3 or genbox.
> c
> c Arguments:
> c
> c uz(nx1,ny1,nlxy): extracted z-slice data
> c u (nx1,ny1,nz1,nelt): input data
> c w1(nx1,ny1,nlxy): work array
> c kz: z-plane within element ezi to be
> extracted
> c ezi: elemental z-slab to be extracted
> c nx,ny,nz: dimensions for 3D spectral element input
> c nlxy: global number of elements in x-y plane.
>
> include 'SIZE'
> include 'GEOM'
> include 'PARALLEL'
> include 'WZ'
>
> real uz(nx,ny,nlxy),u (nx,ny,nz,nelv),w1(nx,ny,nlxy)
> integer e,eg,ex,ey,ez,ezi
> real dy2
>
> nxy = nx*ny
> mxy = nxy*nlxy
>
> call rzero(uz,mxy) ! zero out the entire plane
>
> do e=1,nelt
> eg = lglel(e)
> call get_exyz(ex,ey,ez,eg,nlxy,1,1)
>
> if (ez.eq.ezi)
> $ call copy(uz(1,1,ex),u(1,1,kz,e),nxy) ! u(zlice) --> uz()
>
> enddo
>
> call gop(uz,w1,'+ ',mxy) ! Collect partial contributions from
> all procs
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine z_slice (ua,u,w1,w2)
> c
> c Extract a z slice of quantity u() - assumes global tens.prod.
> c
> include 'SIZE'
> include 'GEOM'
> include 'PARALLEL'
> include 'WZ'
> include 'ZPER'
> c
> real ua(nx1,ny1,nelx,nely),u (nx1,ny1,nz1,nelv)
> $ ,w1(nx1,ny1,nelx,nely),w2(nx1,ny1,nelx,nely)
> integer e,eg,ex,ey,ez
> real dy2
> c
> mxy = nelx*nely*nx1*ny1
> call rzero(ua,mxy)
> c
> do e=1,nelt
> c
> eg = lglel(e)
> call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
>
> k = 1
> if (ez.eq.1) then
> do j=1,ny1
> do i=1,nx1
> ua(i,j,ex,ey) = u(i,j,k,e)
> enddo
> enddo
> endif
> enddo
>
> call gop(ua,w2,'+ ',mxy)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine y_average(ua,u,w1,w2)
> c
> c Compute the y average of quantity u() - assumes global tens.prod.
> c
> include 'SIZE'
> include 'GEOM'
> include 'PARALLEL'
> include 'WZ'
> include 'ZPER'
> c
> real ua(nx1,nz1,nelx,nelz),u (nx1,ny1,nz1,nelv)
> $ ,w1(nx1,nz1,nelx,nelz),w2(nx1,nz1,nelx,nelz)
> integer e,eg,ex,ey,ez
> real dy2
> c
> mxz = nelx*nelz*nx1*nz1
> call rzero(ua,mxz)
> call rzero(w1,mxz)
> c
> do e=1,nelt
> c
> eg = lglel(e)
> call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
> c
> do k=1,nz1
> do i=1,nx1
> c dy2 = .5*( ym1(i,ny1,k,e) - ym1(i,1,k,e) )
> dy2 = 1.0 ! Assuming uniform in "y" direction
> do j=1,ny1
> ua(i,k,ex,ez) = ua(i,k,ex,ez)+dy2*wym1(j)*u(i,j,k,e)
> w1(i,k,ex,ez) = w1(i,k,ex,ez)+dy2*wym1(j) ! redundant
> but clear
> enddo
> enddo
> enddo
> enddo
> c
> call gop(ua,w2,'+ ',mxz)
> call gop(w1,w2,'+ ',mxz)
> c
> do i=1,mxz
> ua(i,1,1,1) = ua(i,1,1,1) / w1(i,1,1,1) ! Normalize
> enddo
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine y_avg_buff(ux,uy,uz,c2,name,icount)
> c
> c Compute the y average of quantity u() - assumes global tens.prod.
> c
> include 'SIZE'
> include 'TOTAL'
> include 'ZPER'
>
> real ux(1),uy(1),uz(1)
> character*2 c2,name
>
> parameter (lyavg = lx1*lz1*lelx*lelz)
> common /scravg/ u (lyavg)
> $ , v (lyavg)
> $ , w (lyavg)
> $ , w1(lyavg)
> $ , w2(lyavg)
>
> call y_average(u,ux,w1,w2)
> call y_average(v,uy,w1,w2)
> call y_average(w,uz,w1,w2)
>
> call buff_2d_out(u,v,w,nx1,nz1,nelx,nelz,c2,name,icount)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine z_avg_buff(ux,uy,uz,c2,name,icount)
> c
> c Compute the z average of quantity u() - assumes global tens.prod.
> c
> include 'SIZE'
> include 'TOTAL'
> include 'ZPER'
> c
> real ux(1),uy(1),uz(1)
> character*2 c2,name
> c
> parameter (lyavg = lx1*ly1*lelx*lely)
> common /scravg/ u (lyavg)
> $ , v (lyavg)
> $ , w (lyavg)
> $ , w1(lyavg)
> $ , w2(lyavg)
> c
> call z_average(u,ux,w1,w2)
> call z_average(v,uy,w1,w2)
> call z_average(w,uz,w1,w2)
>
> call buff_2d_out(u,v,w,nx1,ny1,nelx,nely,c2,name,icount)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine y_ins_buff(ux,uy,uz,c2,name,icount)
> c
> c Compute the z average of quantity u() - assumes global tens.prod.
> c
> include 'SIZE'
> include 'TOTAL'
> include 'ZPER'
> c
> real ux(1),uy(1),uz(1)
> character*2 c2,name
> c
> parameter (lyavg = lx1*lz1*lelx*lelz)
> common /scravg/ u (lyavg)
> $ , v (lyavg)
> $ , w (lyavg)
> $ , w1(lyavg)
> $ , w2(lyavg)
> c
> call y_slice (u,ux,w1,w2)
> call y_slice (v,uy,w1,w2)
> call y_slice (w,uz,w1,w2)
> c
> call buff_2d_out(u,v,w,nx1,nz1,nelx,nelz,c2,name,icount)
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine z_ins_buff(ux,uy,uz,c2,name,icount)
> c
> c Compute the z average of quantity u() - assumes global tens.prod.
> c
> include 'SIZE'
> include 'TOTAL'
> include 'ZPER'
> c
> real ux(1),uy(1),uz(1)
> character*2 c2,name
> c
> parameter (lyavg = lx1*ly1*lelx*lely)
> common /scravg/ u (lyavg)
> $ , v (lyavg)
> $ , w (lyavg)
> $ , w1(lyavg)
> $ , w2(lyavg)
> c
> call z_slice (u,ux,w1,w2)
> call z_slice (v,uy,w1,w2)
> call z_slice (w,uz,w1,w2)
> c
> call buff_2d_out(u,v,w,nx1,ny1,nelx,nely,c2,name,icount)
> c
> return
> end
> c-----------------------------------------------------------------------
> subroutine buff_2d_out(u,v,w,nx,ny,nex,ney,c2,name,ifld)
> c
> include 'SIZE'
> include 'TOTAL'
>
> real u(1),v(1),w(1)
> character*2 c2,name
> character*4 bname
> save bname
>
> parameter (lyzm = lelx*max(lely,lelz))
> common /scrav2/
> ub(lx1,lz1,lyzm),vb(lx1,lz1,lyzm),wb(lx1,lz1,lyzm)
>
> integer ibfld,icalld,nxf,nyf,nexf,neyf
> save ibfld,icalld,nxf,nyf,nexf,neyf
> data ibfld,icalld,nxf,nyf,nexf,neyf / 6*0 /
>
> c npido = 64 ! 64 files buffered
> npido = 128 ! 64 files buffered
> npido = min(npido,np) ! P files buffered
>
> mpido = np/npido ! stride between processors (e.g., 128/64 =
> 2)
>
> jcalld = mod(icalld,npido) ! call # 0,1,...,63,0,1,...
> if (mod(nid,mpido) .eq. 0) then ! this is a buffering/writing
> proc
>
> jid = nid/mpido
> if (jid.eq.jcalld) then ! save this buffer on this proc
> c write(6,1) nid,jid,istep,icalld,jcalld,c2,name,nex,ney,ifld
> c 1 format(5i7,' buffering: ',2a2,3i7)
> write(bname,4) c2,name
> 4 format(2a2)
> n = nx*ny*nex*ney
> ibfld = ifld
> call copy(ub,u,n)
> call copy(vb,v,n)
> call copy(wb,w,n)
> nxf = nx
> nyf = ny
> nexf = nex
> neyf = ney
> endif
>
> if (jcalld .eq. npido-1) call ! output buffer
> $ outfld2d_p(ub,vb,wb,nxf,nyf,nexf,neyf,bname,ibfld,jid,npido,ir)
>
> endif
> call err_chk(ir,'Error with byte_write, buff_2d_out $')
> icalld = icalld+1
> return
> end
> c-----------------------------------------------------------------------
> subroutine y2d(u,v,w,p,c1,icount)
> c
> c Compute the y average of quantity u() - assumes global tens.prod.
> c
>
> include 'SIZE'
> include 'TOTAL'
> real u(1),v(1),w(1),p(1)
> character*1 c1,c2(2)
>
> common /scrns/ ur(lx1*ly1*lz1*lelv)
> $ , ut(lx1*ly1*lz1*lelv)
> $ , wr(lx1*ly1*lz1*lelv)
> $ , wt(lx1*ly1*lz1*lelv)
> $ , wp(lx1*ly1*lz1*lelv)
> c
> c Convert velocities to poloidal-toroidal
> c
> n = nx1*ny1*nz1*nelv
> do i=1,n
> rr = xm1(i,1,1,1)*xm1(i,1,1,1)+ym1(i,1,1,1)*ym1(i,1,1,1)
> rr = sqrt(rr)
> ct = xm1(i,1,1,1)/rr
> st = ym1(i,1,1,1)/rr
> ur(i) = ct*u(i)+st*v(i)
> ut(i) = ct*v(i)-st*u(i)
> wr(i) = ur(i)**2
> wt(i) = ut(i)**2
> wp(i) = w (i)**2
> enddo
>
> c2(1) = c1
> c2(2) = 'y'
>
> call y_avg_buff(ur,w ,ut,c2,'ub',icount)
> call y_avg_buff(wr,wp,wt,c2,'u2',icount)
>
> do i=1,n
> wr(i) = ur(i)*ut(i)
> wt(i) = ut(i)*w (i)
> wp(i) = w (i)*ur(i)
> enddo
> call y_avg_buff(wr,wt,wp,c2,'uv',icount)
>
> call y_ins_buff(ur,w ,ut,c2,'ui',icount)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine z2d(u,v,w,p,c1,icount)
> c
> c Compute the y average of quantity u() - assumes global tens.prod.
> c
>
> include 'SIZE'
> include 'TOTAL'
> real u(1),v(1),w(1),p(1)
> character*1 c1,c2(2)
>
> common /scrns/ ur(lx1*ly1*lz1*lelv)
> $ , ut(lx1*ly1*lz1*lelv)
> $ , wr(lx1*ly1*lz1*lelv)
> $ , wt(lx1*ly1*lz1*lelv)
> $ , wp(lx1*ly1*lz1*lelv)
> c
> c
> c Convert velocities to poloidal-toroidal
> c
> n = nx1*ny1*nz1*nelv
> do i=1,n
> wr(i) = u (i)**2
> wt(i) = v (i)**2
> wp(i) = w (i)**2
> enddo
>
> c2(1) = c1
> c2(2) = 'z'
>
> call z_avg_buff(u ,v ,w ,c2,'ub',icount)
> call z_avg_buff(wr,wt,wp,c2,'u2',icount)
>
> do i=1,n
> wr(i) = u(i)*v(i)
> wt(i) = v(i)*w(i)
> wp(i) = w(i)*u(i)
> enddo
> call z_avg_buff(wr,wt,wp,c2,'uv',icount)
>
> call z_ins_buff(u ,v ,w ,c2,'ui',icount)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine anal_2d
>
> include 'SIZE'
> include 'TOTAL'
> include 'ZPER'
>
> integer icount
> save icount
>
> if (nelx.gt.lelx .or.
> $ nely.gt.lely .or.
> $ nelz.gt.lelz ) then
> if (nid.eq.0) write(6,1) nelx,nely,nelz,lelx,lely,lelz
> 1 format('anal_2d fail:',6i6)
> return
> endif
>
> if (istep.eq.0) then ! dump four times, just to keep phase
>
> icount = 0
> call z2d(xm1,ym1,zm1,pr,'u',icount)
> if (ifmhd) call z2d(xm1,ym1,zm1,pm,'b',icount)
>
> call y2d(xm1,ym1,zm1,pr,'u',icount)
> if (ifmhd) call y2d(xm1,ym1,zm1,pm,'b',icount)
>
> endif
>
> icount = icount + 1
>
> call z2d(vx,vy,vz,pr,'u',icount)
> if (ifmhd) call z2d(bx,by,bz,pm,'b',icount)
>
> call y2d(vx,vy,vz,pr,'u',icount)
> if (ifmhd) call y2d(bx,by,bz,pm,'b',icount)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine chkit(u,name4,n)
>
> include 'SIZE'
> include 'TOTAL'
>
> character*4 name4
> real u(1)
>
>
> integer icalld
> save icalld
> data icalld /0/
>
> icalld = icalld + 1
>
> u2 = vlsc2(u,u,n)
> umin = vlmin(u,n)
> umax = vlmax(u,n)
> ulst = u(n)
> if (nid.eq.0)
> $write(6,1) nid,icalld,istep,n,umin,umax,ulst,name4,' chkit',nid
> 1 format(4i7,1p3e12.4,a4,a6,i1)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine outmesh
> include 'SIZE'
> include 'TOTAL'
> integer e,eg
>
> common /cmesh/ xt(2**ldim,ldim)
>
> len = wdsize*ndim*(2**ndim)
>
> if (nid.eq.0) open(unit=29,file='rea.new')
>
> do eg=1,nelgt
> mtype = eg
> call nekgsync() ! belt
> jnid = gllnid(eg)
> e = gllel (eg)
> if (jnid.eq.0 .and. nid.eq.0) then
> call get_el(xt,xm1(1,1,1,e),ym1(1,1,1,e),zm1(1,1,1,e))
> call out_el(xt,eg)
> elseif (nid.eq.0) then
> call crecv(mtype,xt,len)
> call out_el(xt,eg)
> elseif (jnid.eq.nid) then
> call get_el(xt,xm1(1,1,1,e),ym1(1,1,1,e),zm1(1,1,1,e))
> call csend(mtype,xt,len,0,0)
> endif
> call nekgsync() ! suspenders
> enddo
>
> if (nid.eq.0) close(29)
> call nekgsync()
> call exitt
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine out_el(xt,e)
> include 'SIZE'
> include 'TOTAL'
>
> real xt(2**ldim,ldim)
> integer e
>
> integer ed(8)
> save ed
> data ed / 1,2,4,3 , 5,6,8,7 /
>
> write(29,1) e
> write(29,2) ((xt(ed(k),j),k=1,4),j=1,ndim)
> write(29,2) ((xt(ed(k),j),k=5,8),j=1,ndim)
>
> 1 format(12x,'ELEMENT',i6,' [ 1 ] GROUP 0')
> 2 format(1p4e18.10)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine get_el(xt,x,y,z)
> include 'SIZE'
> include 'TOTAL'
>
> real xt(2**ldim,ldim)
> real x(nx1,ny1,nz1),y(nx1,ny1,nz1),z(nx1,ny1,nz1)
>
> l = 0
> do k=1,nz1,nz1-1
> do j=1,ny1,ny1-1
> do i=1,nx1,nx1-1
> l = l+1
> xt(l,1) = x(i,j,k)
> xt(l,2) = y(i,j,k)
> xt(l,3) = z(i,j,k)
> enddo
> enddo
> enddo
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine shear_calc_max(strsmx,scale,x0,ifdout,iftout)
> c
> c Compute maximum shear stress on objects
> c
> c Scale is a user-supplied multiplier so that results may be
> c scaled to any convenient non-dimensionalization.
> c
> c
> INCLUDE 'SIZE'
> INCLUDE 'TOTAL'
>
> real strsmx(maxobj),x0(3),w1(0:maxobj)
> logical ifdout,iftout
>
> common /cvflow_r/ flow_rate,base_flow,domain_length,xsec
> $ , scale_vf(3)
>
>
> common /scrns/ sij (lx1*ly1*lz1*6*lelv)
> common /scrcg/ pm1 (lx1,ly1,lz1,lelv)
> common /scrsf/ xm0(lx1,ly1,lz1,lelt)
> $, ym0(lx1,ly1,lz1,lelt)
> $, zm0(lx1,ly1,lz1,lelt)
>
> parameter (lr=lx1*ly1*lz1)
> common /scruz/ ur(lr),us(lr),ut(lr)
> $ , vr(lr),vs(lr),vt(lr)
> $ , wr(lr),ws(lr),wt(lr)
>
>
> n = nx1*ny1*nz1*nelv
>
> call mappr(pm1,pr,xm0,ym0) ! map pressure onto Mesh 1
>
> c Add mean_pressure_gradient.X to p:
>
> if (param(55).ne.0) then
> dpdx_mean = -scale_vf(1)
> dpdy_mean = -scale_vf(2)
> dpdz_mean = -scale_vf(3)
> endif
>
> call add2s2(pm1,xm1,dpdx_mean,n) ! Doesn't work if object is cut
> by
> call add2s2(pm1,ym1,dpdy_mean,n) ! periodicboundary. In this
> case,
> call add2s2(pm1,zm1,dpdz_mean,n) ! set ._mean=0 and compensate
> in
> c
> c Compute sij
> c
> nij = 3
> if (if3d.or.ifaxis) nij=6
> call comp_sij(sij,nij,vx,vy,vz,ur,us,ut,vr,vs,vt,wr,ws,wt)
> c
> c
> c Fill up viscous array w/ default
> c
> if (istep.lt.1) call cfill(vdiff,param(2),n)
> c
> call cadd2(xm0,xm1,-x0(1),n)
> call cadd2(ym0,ym1,-x0(2),n)
> call cadd2(zm0,zm1,-x0(3),n)
> c
> x1min=glmin(xm0(1,1,1,1),n)
> x2min=glmin(ym0(1,1,1,1),n)
> x3min=glmin(zm0(1,1,1,1),n)
> c
> x1max=glmax(xm0(1,1,1,1),n)
> x2max=glmax(ym0(1,1,1,1),n)
> x3max=glmax(zm0(1,1,1,1),n)
> c
> call rzero(strsmx,maxobj)
> c
> c
> nobj = 0
> do ii=1,nhis
> if (hcode(10,ii).EQ.'I') then
> iobj = lochis(1,ii)
> memtot = nmember(iobj)
> nobj = max(iobj,nobj)
> c
> if (hcode(1,ii).ne.' ' .or. hcode(2,ii).ne.' ' .or.
> $ hcode(3,ii).ne.' ' ) then
> ifield = 1
> c
> c Compute max stress for this object
> c
> strsmx(ii) = 0.
> do mem=1,memtot
> ieg = object(iobj,mem,1)
> ifc = object(iobj,mem,2)
> if (gllnid(ieg).eq.nid) then ! this processor has a
> contribution
>
> ie = gllel(ieg)
> call get_strsmax
> $ (strsmxl,xm0,ym0,zm0,sij,pm1,vdiff,ifc,ie)
>
> call cmult(strsmxl,scale,1)
> strsmx(ii)=max(strsmx(ii),strsmxl)
>
> endif
> enddo
> endif
> endif
> enddo
> c
> c Max contributions over all processors
> c
> call gop(strsmx,w1,'M ',maxobj)
>
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine get_strsmax(strsmax,xm0,ym0,zm0,sij,pm1,visc,f,e)
> c
> INCLUDE 'SIZE'
> INCLUDE 'GEOM'
> INCLUDE 'INPUT'
> INCLUDE 'TOPOL'
> INCLUDE 'TSTEP'
> c
> real dgtq(3,4)
> real xm0 (lx1,ly1,lz1,lelt)
> real ym0 (lx1,ly1,lz1,lelt)
> real zm0 (lx1,ly1,lz1,lelt)
> real sij (lx1,ly1,lz1,3*ldim-3,lelv)
> real pm1 (lx1,ly1,lz1,lelv)
> real visc(lx1,ly1,lz1,lelv)
>
> integer f,e,pf
> real n1,n2,n3
>
> call dsset(nx1,ny1,nz1) ! set up counters
> pf = eface1(f) ! convert from preproc. notation
> js1 = skpdat(1,pf)
> jf1 = skpdat(2,pf)
> jskip1 = skpdat(3,pf)
> js2 = skpdat(4,pf)
> jf2 = skpdat(5,pf)
> jskip2 = skpdat(6,pf)
>
> if (if3d.or.ifaxis) then
> i = 0
> strsmax = 0
> do j2=js2,jf2,jskip2
> do j1=js1,jf1,jskip1
> i = i+1
> n1 = unx(i,1,f,e)
> n2 = uny(i,1,f,e)
> n3 = unz(i,1,f,e)
> c
> v = visc(j1,j2,1,e)
> c
> s11 = sij(j1,j2,1,1,e)
> s21 = sij(j1,j2,1,4,e)
> s31 = sij(j1,j2,1,6,e)
> c
> s12 = sij(j1,j2,1,4,e)
> s22 = sij(j1,j2,1,2,e)
> s32 = sij(j1,j2,1,5,e)
>
> s13 = sij(j1,j2,1,6,e)
> s23 = sij(j1,j2,1,5,e)
> s33 = sij(j1,j2,1,3,e)
>
> stress1 = -v*(s11*n1 + s12*n2 + s13*n3) ! viscous drag
> stress2 = -v*(s21*n1 + s22*n2 + s23*n3)
> stress3 = -v*(s31*n1 + s32*n2 + s33*n3)
>
> strsnrm = stress1*stress1+stress2*stress2+stress3*stress3
> strsmax = max(strsmax,strsnrm)
>
> enddo
> enddo
>
> else ! 2D
>
> i = 0
> strsmax = 0
> do j2=js2,jf2,jskip2
> do j1=js1,jf1,jskip1
> i = i+1
> n1 = unx(i,1,f,e)*area(i,1,f,e)
> n2 = uny(i,1,f,e)*area(i,1,f,e)
> a = a + area(i,1,f,e)
> v = visc(j1,j2,1,e)
>
> s11 = sij(j1,j2,1,1,e)
> s12 = sij(j1,j2,1,3,e)
> s21 = sij(j1,j2,1,3,e)
> s22 = sij(j1,j2,1,2,e)
>
> stress1 = -v*(s11*n1 + s12*n2) ! viscous drag
> stress2 = -v*(s21*n1 + s22*n2)
>
> strsnrm = stress1*stress1+stress2*stress2
> strsmax = max(strsmax,strsnrm)
>
> enddo
> enddo
>
> endif
>
> if (strsmax.gt.0) strsmax = sqrt(strsmax)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine fix_geom ! fix up geometry irregularities
>
> include 'SIZE'
> include 'TOTAL'
>
> parameter (lt = lx1*ly1*lz1)
> common /scrns/ xb(lt,lelt),yb(lt,lelt),zb(lt,lelt)
> common /scruz/ tmsk(lt,lelt),tmlt(lt,lelt),w1(lt),w2(lt)
>
> integer e,f
> character*3 cb
>
> n = nx1*ny1*nz1*nelt
> nxyz = nx1*ny1*nz1
> nfaces = 2*ndim
> ifield = 1 ! velocity field
> if (ifheat) ifield = 2 ! temperature field
>
>
> call rone (tmlt,n)
> call dssum (tmlt,nx1,ny1,nz1) ! denominator
>
> call rone (tmsk,n)
> do e=1,nelfld(ifield) ! fill mask where bc is periodic
> do f=1,nfaces ! so we don't translate periodic bcs
> (z only)
> cb =cbc(f,e,ifield)
> if (cb.eq.'P ') call facev (tmsk,e,f,0.0,nx1,ny1,nz1)
> enddo
> enddo
>
> do kpass = 1,ndim+1 ! This doesn't work for 2D, yet.
> ! Extra pass is just to test convergence
>
> c call opcopy (xb,yb,zb,xm1,ym1,zm1) ! Must use WHOLE field,
> c call opdssum(xb,yb,zb) ! not just fluid domain.
> call copy (xb,xm1,n)
> call copy (yb,ym1,n)
> call copy (zb,zm1,n)
> call dssum (xb,nx1,ny1,nz1)
> call dssum (yb,nx1,ny1,nz1)
> call dssum (zb,nx1,ny1,nz1)
>
> xm = 0.
> ym = 0.
> zm = 0.
>
> do e=1,nelfld(ifield)
> do i=1,nxyz ! compute averages of
> geometry
> s = 1./tmlt(i,e)
> xb(i,e) = s*xb(i,e)
> yb(i,e) = s*yb(i,e)
> zb(i,e) = s*zb(i,e)
>
> xb(i,e) = xb(i,e) - xm1(i,1,1,e) ! local displacements
> yb(i,e) = yb(i,e) - ym1(i,1,1,e)
> zb(i,e) = zb(i,e) - zm1(i,1,1,e)
> xb(i,e) = xb(i,e)*tmsk(i,e)
> yb(i,e) = yb(i,e)*tmsk(i,e)
> zb(i,e) = zb(i,e)*tmsk(i,e)
>
> xm = max(xm,abs(xb(i,e)))
> ym = max(ym,abs(yb(i,e)))
> zm = max(zm,abs(zb(i,e)))
> enddo
>
> if (kpass.le.ndim) then
> call gh_face_extend(xb(1,e),zgm1,nx1,kpass,w1,w2)
> call gh_face_extend(yb(1,e),zgm1,nx1,kpass,w1,w2)
> call gh_face_extend(zb(1,e),zgm1,nx1,kpass,w1,w2)
> endif
>
> enddo
>
> if (kpass.le.ndim) then
> call add2(xm1,xb,n)
> call add2(ym1,yb,n)
> call add2(zm1,zb,n)
> endif
>
> xx = glamax(xb,n)
> yx = glamax(yb,n)
> zx = glamax(zb,n)
>
> xm = glmax(xm,1)
> ym = glmax(ym,1)
> zm = glmax(zm,1)
>
> if (nid.eq.0) write(6,1) xm,ym,zm,xx,yx,zx,kpass
> 1 format(1p6e12.4,' xyz repair',i2)
>
> enddo
>
> param(59) = 1. ! ifdef = .true.
> call geom_reset(1) ! reset metrics, etc.
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine gh_face_extend(x,zg,n,gh_type,e,v)
> include 'SIZE'
>
> real x(1),zg(1),e(1),v(1)
> integer gh_type
>
> if (ndim.eq.2) then
> call gh_face_extend_2d(x,zg,n,gh_type,e,v)
> else
> call gh_face_extend_3d(x,zg,n,gh_type,e,v)
> endif
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine gh_face_extend_2d(x,zg,n,gh_type,e,v)
> c
> c Extend 2D faces into interior via gordon hall
> c
> c gh_type: 1 - vertex only
> c 2 - vertex and faces
> c
> c
> real x(n,n)
> real zg(n)
> real e(n,n)
> real v(n,n)
> integer gh_type
> c
> c Build vertex interpolant
> c
> ntot=n*n
> call rzero(v,ntot)
> do jj=1,n,n-1
> do ii=1,n,n-1
> do j=1,n
> do i=1,n
> si = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
> sj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
> v(i,j) = v(i,j) + si*sj*x(ii,jj)
> enddo
> enddo
> enddo
> enddo
> if (gh_type.eq.1) then
> call copy(x,v,ntot)
> return
> endif
>
>
> c Extend 4 edges
> call rzero(e,ntot)
> c
> c x-edges
> c
> do jj=1,n,n-1
> do j=1,n
> do i=1,n
> hj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
> e(i,j) = e(i,j) + hj*(x(i,jj)-v(i,jj))
> enddo
> enddo
> enddo
> c
> c y-edges
> c
> do ii=1,n,n-1
> do j=1,n
> do i=1,n
> hi = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
> e(i,j) = e(i,j) + hi*(x(ii,j)-v(ii,j))
> enddo
> enddo
> enddo
>
> call add3(x,e,v,ntot)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine gh_face_extend_3d(x,zg,n,gh_type,e,v)
> c
> c Extend faces into interior via gordon hall
> c
> c gh_type: 1 - vertex only
> c 2 - vertex and edges
> c 3 - vertex, edges, and faces
> c
> c
> real x(n,n,n)
> real zg(n)
> real e(n,n,n)
> real v(n,n,n)
> integer gh_type
> c
> c Build vertex interpolant
> c
> ntot=n*n*n
> call rzero(v,ntot)
> do kk=1,n,n-1
> do jj=1,n,n-1
> do ii=1,n,n-1
> do k=1,n
> do j=1,n
> do i=1,n
> si = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
> sj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
> sk = 0.5*((n-kk)*(1-zg(k))+(kk-1)*(1+zg(k)))/(n-1)
> v(i,j,k) = v(i,j,k) + si*sj*sk*x(ii,jj,kk)
> enddo
> enddo
> enddo
> enddo
> enddo
> enddo
> if (gh_type.eq.1) then
> call copy(x,v,ntot)
> return
> endif
> c
> c
> c Extend 12 edges
> call rzero(e,ntot)
> c
> c x-edges
> c
> do kk=1,n,n-1
> do jj=1,n,n-1
> do k=1,n
> do j=1,n
> do i=1,n
> hj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
> hk = 0.5*((n-kk)*(1-zg(k))+(kk-1)*(1+zg(k)))/(n-1)
> e(i,j,k) = e(i,j,k) + hj*hk*(x(i,jj,kk)-v(i,jj,kk))
> enddo
> enddo
> enddo
> enddo
> enddo
> c
> c y-edges
> c
> do kk=1,n,n-1
> do ii=1,n,n-1
> do k=1,n
> do j=1,n
> do i=1,n
> hi = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
> hk = 0.5*((n-kk)*(1-zg(k))+(kk-1)*(1+zg(k)))/(n-1)
> e(i,j,k) = e(i,j,k) + hi*hk*(x(ii,j,kk)-v(ii,j,kk))
> enddo
> enddo
> enddo
> enddo
> enddo
> c
> c z-edges
> c
> do jj=1,n,n-1
> do ii=1,n,n-1
> do k=1,n
> do j=1,n
> do i=1,n
> hi = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
> hj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
> e(i,j,k) = e(i,j,k) + hi*hj*(x(ii,jj,k)-v(ii,jj,k))
> enddo
> enddo
> enddo
> enddo
> enddo
> c
> call add2(e,v,ntot)
> c
> if (gh_type.eq.2) then
> call copy(x,e,ntot)
> return
> endif
> c
> c Extend faces
> c
> call rzero(v,ntot)
> c
> c x-edges
> c
> do ii=1,n,n-1
> do k=1,n
> do j=1,n
> do i=1,n
> hi = 0.5*((n-ii)*(1-zg(i))+(ii-1)*(1+zg(i)))/(n-1)
> v(i,j,k) = v(i,j,k) + hi*(x(ii,j,k)-e(ii,j,k))
> enddo
> enddo
> enddo
> enddo
> c
> c y-edges
> c
> do jj=1,n,n-1
> do k=1,n
> do j=1,n
> do i=1,n
> hj = 0.5*((n-jj)*(1-zg(j))+(jj-1)*(1+zg(j)))/(n-1)
> v(i,j,k) = v(i,j,k) + hj*(x(i,jj,k)-e(i,jj,k))
> enddo
> enddo
> enddo
> enddo
> c
> c z-edges
> c
> do kk=1,n,n-1
> do k=1,n
> do j=1,n
> do i=1,n
> hk = 0.5*((n-kk)*(1-zg(k))+(kk-1)*(1+zg(k)))/(n-1)
> v(i,j,k) = v(i,j,k) + hk*(x(i,j,kk)-e(i,j,kk))
> enddo
> enddo
> enddo
> enddo
> c
> call add2(v,e,ntot)
> call copy(x,v,ntot)
>
> return
> end
> c-----------------------------------------------------------------------
> function ran1(idum)
> c
> integer idum,ia,im,iq,ir,ntab,ndiv
> real ran1,am,eps,rnmx
> c
> parameter (ia=16807,im=2147483647,am=1./im,iq=127773,ir=2836)
> parameter (ntab=32,ndiv=1+(im-1)/ntab,eps=1.2e-7,rnmx=1.-eps)
> c
> c Numerical Rec. in Fortran, 2nd eD. P. 271
> c
> integer j,k
> integer iv(ntab),iy
> save iv,iy
> data iv,iy /ntab*0,0/
> c
> if (idum.le.0.or.iy.eq.0) then
> idum=max(-idum,1)
> do j=ntab+8,1,-1
> k = idum/iq
> idum = ia*(idum-k*iq)-ir*k
> if(idum.lt.0) idum = idum+im
> if (j.le.ntab) iv(j) = idum
> enddo
> iy = iv(1)
> endif
> k = idum/iq
> idum = ia*(idum-k*iq)-ir*k
> if(idum.lt.0) idum = idum+im
> j = 1+iy/ndiv
> iy = iv(j)
> iv(j) = idum
> ran1 = min(am*iy,rnmx)
> c ran1 = cos(ran1*1.e8)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine rand_fld_h1(x)
>
> include 'SIZE'
> real x(1)
>
> n=nx1*ny1*nz1*nelt
> id = n
> do i=1,n
> x(i) = ran1(id)
> enddo
> call dsavg(x)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine rescale_x (x,x0,x1)
> include 'SIZE'
> real x(1)
>
> n = nx1*ny1*nz1*nelt
> xmin = glmin(x,n)
> xmax = glmax(x,n)
>
> if (xmax.le.xmin) return
>
> scale = (x1-x0)/(xmax-xmin)
> do i=1,n
> x(i) = x0 + scale*(x(i)-xmin)
> enddo
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine z_distribute(u)
> c
> c Compute the z average of quantity u() and redistribute
> c
> c Assumes you have nelx*nely elements, in the same order,
> c within each z plane
> c
> c
> include 'SIZE'
> include 'TOTAL'
> include 'ZPER'
>
> real ux(1),uy(1),uz(1)
> character*2 c2,name
>
> parameter (lyavg = lx1*ly1*lelx*lely)
> common /scravg/ ua(lyavg)
> $ , w1(lyavg)
> $ , w2(lyavg)
>
> call z_average (ua,u,w1,w2)
> call z_average_transpose(u,ua) ! distribute ua to each z-plane
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine z_average(ua,u,w1,w2)
> c
> c Compute the z average of quantity u() - assumes global tens.prod.
> c
> include 'SIZE'
> include 'GEOM'
> include 'PARALLEL'
> include 'WZ'
> include 'ZPER'
>
> real ua(nx1,ny1,nelx,nely),u (nx1,ny1,nz1,nelv)
> $ ,w1(nx1,ny1,nelx,nely),w2(nx1,ny1,nelx,nely)
> integer e,eg,ex,ey,ez
> real dy2
>
> nelxy = nelx*nely
> if (nelxy.gt.lelx*lely) call exitti
> $ ('ABORT IN z_average. Increase lelx*lely in SIZE:$',nelxy)
>
> mxy = nelx*nely*nx1*ny1
> call rzero(ua,mxy)
> call rzero(w1,mxy)
>
> do e=1,nelt
>
> eg = lglel(e)
> call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
>
> do j=1,ny1
> do i=1,nx1
> dz2 = 1.0 ! Assuming uniform in "z" direction
> do k=1,nz1
> ua(i,j,ex,ey) = ua(i,j,ex,ey)+dz2*wzm1(k)*u(i,j,k,e)
> w1(i,j,ex,ey) = w1(i,j,ex,ey)+dz2*wzm1(k) ! redundant
> but clear
> enddo
> enddo
> enddo
> enddo
>
> call gop(ua,w2,'+ ',mxy)
> call gop(w1,w2,'+ ',mxy)
>
> do i=1,mxy
> ua(i,1,1,1) = ua(i,1,1,1) / w1(i,1,1,1) ! Normalize
> enddo
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine z_average_transpose(u,ua) ! distribute ua to each
> z-plane
>
> include 'SIZE'
> include 'GEOM'
> include 'PARALLEL'
> include 'WZ'
> include 'ZPER'
>
> real u(nx1,ny1,nz1,nelv),ua(nx1,ny1,nelx,nely)
>
> integer e,eg,ex,ey,ez
>
>
> do e=1,nelt
>
> eg = lglel(e)
> call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
>
> do j=1,ny1
> do i=1,nx1
> do k=1,nz1
> u(i,j,k,e) = ua(i,j,ex,ey)
> enddo
> enddo
> enddo
> enddo
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine no_z_profile(u)
>
> c Subtract the z_profile from u for a tensor-product array of
> elements
>
> c Assumes you have nelx*nely*nelz elements, in the same order,
> c and that lelx,lely,lelz are defined to be >= nelx,nely,nelz
>
>
> include 'SIZE'
> include 'TOTAL'
> include 'ZPER' ! nelx,nely,nelz
>
> real u(1)
>
> parameter (lyavg = ly1*lely)
> common /scravg/ ua(lyavg)
> $ , w1(lyavg)
> $ , w2(lyavg)
> common /scrmg/ ub(lx1*ly1*lz1*lelt)
>
> call z_profile (ua,u,w1,w2)
> call z_profile_transpose(ub,ua) ! distribute ua to each z-plane
>
> n = nx1*ny1*nz1*nelv
> call sub2(u,ub,n)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine z_profile(ua,u,w1,w2)
> c
> c Compute the z profile of quantity u() - assumes global tens.prod.
> c
> include 'SIZE'
> include 'GEOM'
> include 'PARALLEL'
> include 'WZ'
> include 'ZPER'
>
> real ua(lz1,lelz),u (nx1,ny1,nz1,nelv)
> $ ,w1(lz1,lelz),w2(lz1,lelz)
> integer e,eg,ex,ey,ez
>
> mz = nz1*nelz
> call rzero(ua,mz)
> call rzero(w1,mz)
>
> do e=1,nelt
>
> eg = lglel(e)
> call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
>
> do k=1,nz1
> do i=1,nx1*ny1
> ua(k,ez) = ua(k,ez) + area(i,1,5,e)*u(i,1,k,e)
> w1(k,ez) = w1(k,ez) + area(i,1,5,e)
> enddo
> enddo
>
> enddo
>
> call gop(ua,w2,'+ ',mz)
> call gop(w1,w2,'+ ',mz)
>
> do i=1,mz
> ua(i,1) = ua(i,1) / w1(i,1) ! Normalize
> enddo
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine z_profile_transpose(u,ua) ! distribute ua to each
> z-plane
>
> include 'SIZE'
> include 'PARALLEL'
> include 'ZPER'
>
> real u(nx1,ny1,nz1,nelv),ua(lz1,lelz)
> integer e,eg,ex,ey,ez
>
> do e=1,nelt
>
> eg = lglel(e)
> call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
>
> do k=1,nz1
> do i=1,nx1*ny1
> u(i,1,k,e) = ua(k,ez)
> enddo
> enddo
>
> enddo
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine no_y_profile(u)
>
> c Subtract the y_profile from u for a tensor-product array of
> elements
>
> c Assumes you have nelx*nely*nelz elements, in the same order,
> c and that lelx,lely,lelz are defined to be >= nelx,nely,nelz
>
>
> include 'SIZE'
> include 'TOTAL'
> include 'ZPER' ! nelx,nely,nelz
>
> real u(1)
>
> parameter (lyavg = ly1*lely)
> common /scravg/ ua(lyavg)
> $ , w1(lyavg)
> $ , w2(lyavg)
> common /scrmg/ ub(lx1*ly1*lz1*lelt)
>
> call y_profile (ua,u,w1,w2)
> call y_profile_transpose(ub,ua) ! distribute ua to each y-plane
>
> n = nx1*ny1*nz1*nelv
> call sub2(u,ub,n)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine y_profile(ua,u,w1,w2)
> c
> c Compute the z profile of quantity u() - assumes global tens.prod.
> c
> include 'SIZE'
> include 'GEOM'
> include 'PARALLEL'
> include 'WZ'
> include 'ZPER'
>
> real ua(lz1,lelz),u (nx1,ny1,nz1,nelv)
> $ ,w1(lz1,lelz),w2(lz1,lelz)
> integer e,eg,ex,ey,ez
>
> my = ny1*nely
> call rzero(ua,my)
> call rzero(w1,my)
>
> do e=1,nelt
>
> eg = lglel(e)
> call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
>
> do k=1,nz1
> do j=1,ny1
> do i=1,nx1
> ua(j,ey) = ua(j,ey) + area(i,k,1,e)*u(i,j,k,e)
> w1(j,ey) = w1(j,ey) + area(i,k,1,e)
> enddo
> enddo
> enddo
>
> enddo
>
> call gop(ua,w2,'+ ',my)
> call gop(w1,w2,'+ ',my)
>
> do i=1,my
> ua(i,1) = ua(i,1) / w1(i,1) ! Normalize
> enddo
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine y_profile_transpose(u,ua) ! distribute ua to each
> z-plane
>
> include 'SIZE'
> include 'PARALLEL'
> include 'ZPER'
>
> real u(nx1,ny1,nz1,nelv),ua(lz1,lelz)
> integer e,eg,ex,ey,ez
>
> do e=1,nelt
>
> eg = lglel(e)
> call get_exyz(ex,ey,ez,eg,nelx,nely,nelz)
>
> do k=1,nz1
> do j=1,ny1
> do i=1,nx1
> u(i,j,k,e) = ua(j,ey)
> enddo
> enddo
> enddo
>
> enddo
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine build_filter(f,diag,nx)
> include 'SIZE'
>
> real f(nx,nx),diag(nx),zpts(nx)
>
> parameter (lm=4*lx1) ! Totally arbitrary
> parameter (lm2=lm*lm)
>
> common /cfilt1/ phi,pht,ft,rmult,Lj,gpts,indr,indc,ipiv
> real phi(lm2),pht(lm2),ft(lm2),rmult(lm),Lj(lm),gpts(lm)
> integer indr(lm),indc(lm),ipiv(lm)
>
> integer nxl
> save nxl
> data nxl / -9 /
>
> if (nx.gt.lm) call exitti('ABORT in build_filter:$',nx)
>
> if (nx.ne.nxl) then
>
> nxl = nx
>
> call zwgll (gpts,f,nx) ! Get nx GLL points
>
> kj = 0
> n = nx-1
> do j=1,nx
> z = gpts(j)
> call legendre_poly(Lj,z,n)
> kj = kj+1
> pht(kj) = Lj(1)
> kj = kj+1
> pht(kj) = Lj(2)
> do k=3,nx
> kj = kj+1
> pht(kj) = Lj(k)-Lj(k-2)
> enddo
> enddo
>
> call transpose (phi,nx,pht,nx)
> call copy (pht,phi,nx*nx)
> call gaujordf (pht,nx,nx,indr,indc,ipiv,ierr,rmult)
>
> endif ! End of save section
>
> ij=0
> do j=1,nx
> do i=1,nx
> ij = ij+1
> ft(ij) = diag(i)*pht(ij)
> enddo
> enddo
> ! -1
> call mxm (phi,nx,ft,nx,f,nx) ! V D V
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine g_filter(u,diag,ifld)
> c
> c Generalized filter: F(u) with F = J^T D J, where D=diag(diag)
> c
> include 'SIZE'
> include 'TOTAL'
>
> real u(1),diag(1)
>
> parameter (lxx=lx1*lx1,lxyz=lx1*ly1*lz1)
> common /ctmp0/ f(lxx),wk1(lxyz),wk2(lxyz),wk3(lxyz)
>
> ifldt = ifield
> ifield = ifld
>
> call build_filter(f,diag,nx1)
> call filterq(u,f,nx1,nz1,wk1,wk2,wk3,if3d,umax)
>
> ifield = ifldt
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine cut_off_filter(u,mx,ifld) ! mx=max saved mode
> c
> c Generalized filter: F(u) with F = J^T D J, where D=diag(diag)
> c
> include 'SIZE'
> include 'TOTAL'
>
> real u(1)
>
> parameter (lxx=lx1*lx1,lxyz=lx1*ly1*lz1)
> common /ctmp0/ f(lxx),wk1(lxyz),wk2(lxyz),wk3(lxyz),diag(lx1)
>
> ifldt = ifield
> ifield = ifld
>
> call rone(diag,nx1)
> do i=mx+1,nx1
> diag(i)=0.
> enddo
>
> call build_filter(f,diag,nx1)
> call filterq(u,f,nx1,nz1,wk1,wk2,wk3,if3d,umax)
>
> ifield = ifldt
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine filter_d2(v,nx,nz,wgt,ifd4)
>
> include 'SIZE'
> include 'TOTAL'
>
> parameter (lt=lx1*ly1*lz1)
> real v(lt,nelt)
> logical ifd4
>
> common /ctmp1/ w(lt,lelt),ur(lt),us(lt),ut(lt),w1(2*lt)
>
> integer e
>
> n = nx1*ny1*nz1*nelt
> nn = nx-1
> nel = nelfld(ifield)
>
> bmax = glamax(v,n)
>
> if (if3d) then
> do e=1,nel
> call local_grad3(ur,us,ut,v(1,e),nn,1,dxm1,dxtm1)
> do i=1,lt
> ur(i) = ur(i)*w3m1(i,1,1)
> us(i) = us(i)*w3m1(i,1,1)
> ut(i) = ut(i)*w3m1(i,1,1)
> enddo
> call local_grad3_t(w(1,e),ur,us,ut,nn,1,dxm1,dxtm1,w1)
> enddo
> call dsavg(w) !NOTE STILL NEED BC TREATMENT !
>
> if (ifd4) then
> wght = 20./(nx1**4)
> do e=1,nel
> do i=1,lt
> w(i,e) = wght*w(i,e)/w3m1(i,1,1)
> enddo
> call local_grad3(ur,us,ut,w(1,e),nn,1,dxm1,dxtm1)
> do i=1,lt
> ur(i) = ur(i)*w3m1(i,1,1)
> us(i) = us(i)*w3m1(i,1,1)
> ut(i) = ut(i)*w3m1(i,1,1)
> enddo
> call local_grad3_t(w(1,e),ur,us,ut,nn,1,dxm1,dxtm1,w1)
> enddo
> call dsavg(w) !NOTE STILL NEED BC TREATMENT !
> endif
>
> wght = wgt/(nx1**4)
> do e=1,nel
> do i=1,lt
> v(i,e) = v(i,e) - wght*w(i,e)/w3m1(i,1,1)
> enddo
> enddo
>
> else ! 2D
>
> do e=1,nel
> call local_grad2(ur,us,v(1,e),nn,1,dxm1,dxtm1)
> do i=1,lt
> ur(i) = ur(i)*w3m1(i,1,1)
> us(i) = us(i)*w3m1(i,1,1)
> enddo
> call local_grad2_t(w(1,e),ur,us,nn,1,dxm1,dxtm1,w1)
> enddo
> call dsavg(w) !NOTE STILL NEED BC TREATMENT !
>
> if (ifd4) then
> wght = 200./(nx1**4)
> do e=1,nel
> do i=1,lt
> w(i,e) = wght*w(i,e)/w3m1(i,1,1)
> enddo
> call local_grad2(ur,us,w(1,e),nn,1,dxm1,dxtm1)
> do i=1,lt
> ur(i) = ur(i)*w3m1(i,1,1)
> us(i) = us(i)*w3m1(i,1,1)
> enddo
> call local_grad2_t(w(1,e),ur,us,nn,1,dxm1,dxtm1,w1)
> enddo
> call dsavg(w) !NOTE STILL NEED BC TREATMENT !
> endif
>
> wght = wgt/(nx1**4)
> do e=1,nel
> do i=1,lt
> v(i,e) = v(i,e) - wght*w(i,e)/w3m1(i,1,1)
> enddo
> enddo
>
> endif
>
> vmax = glamax(v,n)
> if (nid.eq.0) write(6,1) istep,time,vmax,bmax,' filter max'
> 1 format(i9,1p3e12.4,a11)
>
> return
> end
> c-------------------------------------------------------------------------
> function dist3d(a,b,c,x,y,z)
>
> d = (a-x)**2 + (b-y)**2 + (c-z)**2
>
> dist3d = 0.
> if (d.gt.0) dist3d = sqrt(d)
>
> return
> end
> c-----------------------------------------------------------------------
> function dist2d(a,b,x,y)
>
> d = (a-x)**2 + (b-y)**2
>
> dist2d = 0.
> if (d.gt.0) dist2d = sqrt(d)
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine domain_size(xmin,xmax,ymin,ymax,zmin,zmax)
>
> include 'SIZE'
> include 'TOTAL'
>
> n = nx1*ny1*nz1*nelt
>
> xmin = glmin(xm1,n)
> xmax = glmax(xm1,n)
> ymin = glmin(ym1,n)
> ymax = glmax(ym1,n)
> if (if3d) then
> zmin = glmin(zm1,n)
> zmax = glmax(zm1,n)
> else
> zmin = 0.
> zmax = 0.
> endif
>
> return
> end
> c-----------------------------------------------------------------------
> subroutine cheap_dist(d,ifld,b)
>
> c Finds a pseudo-distance function.
>
> c INPUT: ifld - field type for which distance function is to be
> found.
> c ifld = 1 for velocity
> c ifld = 2 for temperature, etc.
>
> c OUTPUT: d = "path" distance to nearest wall
>
> c This approach has a significant advantage that it works for
> c periodict boundary conditions, whereas most other approaches
> c will not.
>
> include 'SIZE'
> include 'GEOM' ! Coordinates
> include 'INPUT' ! cbc()
> include 'TSTEP' ! nelfld
> include 'PARALLEL' ! gather-scatter handle for field "ifld"
>
> real d(lx1,ly1,lz1,lelt)
>
> character*3 b ! Boundary condition of interest
>
> integer e,eg,f
>
> nel = nelfld(ifld)
> n = nx1*ny1*nz1*nel
>
> call domain_size(xmin,xmax,ymin,ymax,zmin,zmax)
>
> xmn = min(xmin,ymin)
> xmx = max(xmax,ymax)
> if (if3d) xmn = min(xmn ,zmin)
> if (if3d) xmx = max(xmx ,zmax)
>
> big = 10*(xmx-xmn)
> call cfill(d,big,n)
>
>
> nface = 2*ndim
> do e=1,nel ! Set d=0 on walls
> do f=1,nface
> if (cbc(f,e,ifld).eq.b) call facev(d,e,f,0.,nx1,ny1,nz1)
> enddo
> enddo
>
> do ipass=1,10000
> dmax = 0
> nchange = 0
> do e=1,nel
> do k=1,nz1
> do j=1,ny1
> do i=1,nx1
> i0=max( 1,i-1)
> j0=max( 1,j-1)
> k0=max( 1,k-1)
> i1=min(nx1,i+1)
> j1=min(ny1,j+1)
> k1=min(nz1,k+1)
> do kk=k0,k1
> do jj=j0,j1
> do ii=i0,i1
>
> if (if3d) then
> dtmp = d(ii,jj,kk,e) + dist3d(
> $ xm1(ii,jj,kk,e),ym1(ii,jj,kk,e),zm1(ii,jj,kk,e)
> $ ,xm1(i ,j ,k ,e),ym1(i ,j ,k ,e),zm1(i ,j ,k ,e))
> else
> dtmp = d(ii,jj,kk,e) + dist2d(
> $ xm1(ii,jj,kk,e),ym1(ii,jj,kk,e)
> $ ,xm1(i ,j ,k ,e),ym1(i ,j ,k ,e))
> endif
>
> if (dtmp.lt.d(i,j,k,e)) then
> d(i,j,k,e) = dtmp
> nchange = nchange+1
> dmax = max(dmax,d(i,j,k,e))
> endif
> enddo
> enddo
> enddo
>
> enddo
> enddo
> enddo
> enddo
> call gs_op(gsh_fld(ifld),d,1,3,0) ! min over all elements
> nchange = iglsum(nchange,1)
> dmax = glmax(dmax,1)
> if (nid.eq.0) write(6,1) ipass,nchange,dmax
> 1 format(i9,i12,1pe12.4,' max wall distance 1')
> if (nchange.eq.0) goto 1000
> enddo
> 1000 return
> end
> c-----------------------------------------------------------------------
> subroutine distf(d,ifld,b,dmin,emin,xn,yn,zn)
>
> c Generate a distance function to boundary with bc "b".
> c This approach does not yet work with periodic boundary
> conditions.
>
> c INPUT: ifld - field type for which distance function is to be
> found.
> c ifld = 1 for velocity
> c ifld = 2 for temperature, etc.
>
> c OUTPUT: d = distance to nearest boundary with boundary condition
> "b"
>
> c Work arrays: dmin,emin,xn,yn,zn
>
> include 'SIZE'
> include 'GEOM' ! Coordinates
> include 'INPUT' ! cbc()
> include 'TSTEP' ! nelfld
> include 'PARALLEL' ! gather-scatter handle for field "ifld"
>
> real d(lx1,ly1,lz1,lelt)
> character*3 b
>
> real dmin(lx1,ly1,lz1,lelt),emin(lx1,ly1,lz1,lelt)
> real xn(lx1,ly1,lz1,lelt),yn(lx1,ly1,lz1,lelt)
> real zn(lx1,ly1,lz1,lelt)
>
>
> integer e,eg,f
>
> nel = nelfld(ifld)
> n = nx1*ny1*nz1*nel
>
> call domain_size(xmin,xmax,ymin,ymax,zmin,zmax)
>
> xmn = min(xmin,ymin)
> xmx = max(xmax,ymax)
> if (if3d) xmn = min(xmn ,zmin)
> if (if3d) xmx = max(xmx ,zmax)
>
> big = 10*(xmx-xmn)
> call cfill (d,big,n)
>
> call opcopy(xn,yn,zn,xm1,ym1,zm1)
>
> nface = 2*ndim
> do e=1,nel ! Set d=0 on walls
> do f=1,nface
> if (cbc(f,e,1).eq.b) call facev(d,e,f,0.,nx1,ny1,nz1)
> enddo
> enddo
>
> nxyz = nx1*ny1*nz1
>
> do ipass=1,10000
> dmax = 0
> nchange = 0
> do e=1,nel
> do k=1,nz1
> do j=1,ny1
> do i=1,nx1
> i0=max( 1,i-1)
> j0=max( 1,j-1)
> k0=max( 1,k-1)
> i1=min(nx1,i+1)
> j1=min(ny1,j+1)
> k1=min(nz1,k+1)
> do kk=k0,k1
> do jj=j0,j1
> do ii=i0,i1
>
> dself = d(i,j,k,e)
> dneigh = d(ii,jj,kk,e)
> if (dneigh.lt.dself) then ! check neighbor's nearest
> point
> d2 = (xm1(i,j,k,e)-xn(ii,jj,kk,e))**2
> $ + (ym1(i,j,k,e)-yn(ii,jj,kk,e))**2
> if (if3d) d2 = d2 + (zm1(i,j,k,e)-zn(ii,jj,kk,e))**2
> if (d2.gt.0) d2 = sqrt(d2)
> if (d2.lt.dself) then
> nchange = nchange+1
> d (i,j,k,e) = d2
> xn(i,j,k,e) = xn(ii,jj,kk,e)
> yn(i,j,k,e) = yn(ii,jj,kk,e)
> zn(i,j,k,e) = zn(ii,jj,kk,e)
> dmax = max(dmax,d(i,j,k,e))
> endif
> endif
> enddo
> enddo
> enddo
>
> enddo
> enddo
> enddo
>
> re = lglel(e)
> call cfill(emin(1,1,1,e),re,nxyz)
> call copy (dmin(1,1,1,e),d(1,1,1,e),nxyz)
>
> enddo
> nchange = iglsum(nchange,1)
>
> call gs_op(gsh_fld(ifld),dmin,1,3,0) ! min over all elements
>
>
> nchange2=0
> do e=1,nel
> do i=1,nxyz
> if (dmin(i,1,1,e).ne.d(i,1,1,e)) then
> nchange2 = nchange2+1
> emin(i,1,1,e) = 0 ! Flag
> endif
> enddo
> enddo
> call copy(d,dmin,n) ! Ensure updated distance
> nchange2 = iglsum(nchange2,1)
> nchange = nchange + nchange2
> call gs_op(gsh_fld(ifld),emin,1,4,0) ! max over all elements
>
> do e=1,nel ! Propagate nearest wall points
> do i=1,nxyz
> eg = emin(i,1,1,e)
> if (eg.ne.lglel(e)) then
> xn(i,1,1,e) = 0
> yn(i,1,1,e) = 0
> zn(i,1,1,e) = 0
> endif
> enddo
> enddo
> call gs_op(gsh_fld(ifld),xn,1,1,0) ! Sum over all elements
> to
> call gs_op(gsh_fld(ifld),yn,1,1,0) ! convey nearest point
> call gs_op(gsh_fld(ifld),zn,1,1,0) ! to shared neighbor.
>
> dmax = glmax(dmax,1)
> if (nid.eq.0) write(6,1) ipass,nchange,dmax
> 1 format(i9,i12,1pe12.4,' max wall distance 2')
> if (nchange.eq.0) goto 1000
> enddo
> 1000 continue
>
> c wgt = 0.3
> c call filter_d2(d,nx1,nz1,wgt,.true.)
>
> return
> end
> c-----------------------------------------------------------------------
>
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
More information about the Nek5000-users
mailing list