[Nek5000-users] Extruding 2D .fld file to 3D
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Wed Oct 23 19:52:45 CDT 2013
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-----------------------------------------------------------------------
More information about the Nek5000-users
mailing list