From 1f8402a266d14e99496bd8437d938df0a86f597f Mon Sep 17 00:00:00 2001 From: Dillon Shaver Date: Tue, 11 Dec 2018 11:49:51 -0600 Subject: [PATCH 1/5] use local array for pressure lagging in monitor routine --- RANSChannel/RANSchan.usr | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/RANSChannel/RANSchan.usr b/RANSChannel/RANSchan.usr index 3d75e2c..2cd5f98 100644 --- a/RANSChannel/RANSchan.usr +++ b/RANSChannel/RANSchan.usr @@ -487,6 +487,7 @@ c----------------------------------------------------------------------- real vol,glsum,glmin,glmax,glsc2 real tmp(lx1*ly1*lz1*lelv) + real prlagl(lx2*ly2*lz2*lelv) integer i,n1,n2,nt real uxmin,uxmax,uxave,uymin,uymax,uyave,uzmin,uzmax,uzave @@ -496,6 +497,8 @@ c----------------------------------------------------------------------- real dux,duy,duz,dpr,dth(ldimt) character*15 tname + save prlagl + n1=nx1*ny1*nz1*nelv n2=nx2*ny2*nz2*nelv nt=nx1*ny1*nz1*nelt @@ -503,7 +506,7 @@ c----------------------------------------------------------------------- call get_limits(vx,uxmin,uxmax,uxave,dux,vxlag,rmsux,n1) call get_limits(vy,uymin,uymax,uyave,duy,vylag,rmsuy,n1) if(if3d) call get_limits(vz,uzmin,uzmax,uzave,duz,vzlag,rmsuz,n1) - call get_limits(pr,prmin,prmax,prave,dpr,prlag,rmspr,n2) + call get_limits(pr,prmin,prmax,prave,dpr,prlagl,rmspr,n2) if(ifheat)then do i=1,npscal+1 if(idpss(i).eq.0) call get_limits(t(1,1,1,1,i),thmin(i) !Helmholtz solver @@ -533,7 +536,7 @@ c----------------------------------------------------------------------- write(*,*) endif - call copy(prlag,pr,n2) + call copy(prlagl,pr,n2) 254 format(a15,5a13) 255 format(a15,5es13.4) From b5a39db9c424dda6473da527947bd5592b0908ae Mon Sep 17 00:00:00 2001 From: dshaver-ANL Date: Tue, 2 Mar 2021 22:45:53 -0600 Subject: [PATCH 2/5] update to use implicit k-tau --- RANSChannel/RANSchan.box | 10 +- RANSChannel/RANSchan.ma2 | Bin 616 -> 736 bytes RANSChannel/RANSchan.par | 32 +-- RANSChannel/RANSchan.re2 | Bin 7484 -> 8940 bytes RANSChannel/RANSchan.usr | 555 +++++++++++++++------------------------ RANSChannel/README.md | 2 +- RANSChannel/SIZE | 13 +- 7 files changed, 237 insertions(+), 375 deletions(-) diff --git a/RANSChannel/RANSchan.box b/RANSChannel/RANSchan.box index eaff81e..fe1cd57 100644 --- a/RANSChannel/RANSchan.box +++ b/RANSChannel/RANSchan.box @@ -14,10 +14,10 @@ #======================================================================= # Box --3 -8 nelx,nely,nelz for Box +-3 -10 nelx,nely,nelz for Box 0 1.0 1.0 x0,x1,gain -0.0 1.0 0.55 +0.0 1.0 0.57 P ,P ,SYM,W -P ,P ,SYM,I -P ,P ,SYM,t -P ,P ,SYM,t +P ,P ,I ,f +P ,P ,I ,t +P ,P ,I ,t diff --git a/RANSChannel/RANSchan.ma2 b/RANSChannel/RANSchan.ma2 index 8dbfff728a2cde810afcea047aa28fb41e7abfaf..1a260d4eb16b2c7839af6bf977d969cd0afa784e 100644 GIT binary patch literal 736 zcmZ{hS#H8W6a>u@_OOJV1Q^5x>=^OQRk#HAwqMbo7LV)wGRm^7t}K+bvQmc8=h2s%aCBHY>~~C@nX*%G z=`i%I0%wN~JBSJMr~E2i#SVR=z}cb0c4ET(C=&(O4t=Wh6+3j;UQC#Sk`(sr(Dw?i z9Xf0~X42#wucagfhR%1u*`dRJ#f0Je#|o|;I(H0bhYs6}3B&zzpEbn}oqK|_Lx<(D f!~6JsTsw5Wx1rdf!}4F-;T(@WJ9OU3nH~2hB7YsG literal 616 zcmaiwNe;p=5Cp@%?>k&9@c}17AZ~nxPw-E^!-+0wCN?J`wZ@fH+eyx^7|T#!)umpB z>Xz0ayI)?@{d>0l@ndrSu0GHAV^$PJQ&}n-Wvxt&7se}1JQx?oo!F$aS9Z!+fpOri zg6AC=cVrW1uCV`5@ea&;;CTne9oWP*jZ^R3$aUOp`DbAAS7l*@GSXe-vmlPd{iC_Mu&zbYJi&b)eLPruxE5~Hc$7L&>OoU| z;YXc6&dX+=e$XoHL%TTXzTnYRJ?{SD^537kPn+J)Uplk)myYGPm+pFt*BCP}=2ySQ%jKV!{6{+GSHH%~r^7dcjz*SRq*@t~qwUi}&`a;z+`b3^9}R8-5WU*koNmE{Me zyU+Kfa|J4@<<+n8BFD<|IyZE#Kt;8@`ZZqUSXo}@hRzkJsFqj1#)}*)%j?|GxdIi{ z^6J-kkz-|f&TY;iTxq>w&QExrAI^2oqk0^hzxy1r`8A)`lgkf3ZvPu|$n9|0ztUYZ z&v-6!-vd5W=Sq%OKXK3;uXCecdII&AQS5Wv$LCe!MUIu_b#CZffr@H*^=rJyv9i3* z4V^1cQ7x~2jTbpqmLHVvKHrzl6{x6|SHH%K94pJ~+|aoK71i?U*Laa*WqF+&I#-~g zT3-DcFLJCbuX97^3RG0ft6$?qj+Ny(x5GJtEva|@^Bl7InIp_S`r&!5w_*E-AGiOF zIplUY>|g1wH|CJrhqlk5&JCR_F52k3eU8_ohvzB{TeUF>)g<}lH=8{ z@p8P*4V^1FUi}&`$LrkCxsv16ukmuc&JCR_IbQu5FUNCkhjXM3$Q?Te;Q76At}_<| X{(EVk|L*?}+5DPM>&fMZ@816(ZFl$J delta 960 zcmaFky2naRxy-=8NC5(jOcWH5*cmBB6BYF*`r7l?t`=+LHR8n^AKz z52FBvxHhAvJeqd<=SP;Qy9dU6(r z!o*cln4(q!s0!cPPmbZ2V>Fu_&2PDxhf|DYvJ#)gWGP+;M)S$ieD;&0fGi6T3ss5M gWEn9&xH3a952h5%Lso7i4Kq<&8fGGpMOcX?00Lj7jQ{`u diff --git a/RANSChannel/RANSchan.usr b/RANSChannel/RANSchan.usr index ee0576f..d1e694e 100644 --- a/RANSChannel/RANSchan.usr +++ b/RANSChannel/RANSchan.usr @@ -1,5 +1,4 @@ #include "experimental/rans_komg.f" - c----------------------------------------------------------------------- c nek5000 user-file template c @@ -23,23 +22,27 @@ c----------------------------------------------------------------------- include 'NEKUSE' integer ix,iy,iz,e,eg - real rans_komg_mut,rans_komg_mutsk,rans_komg_mutso + real rans_mut,rans_mutsk,rans_mutso,rans_turbPrandtl real mu_t,Pr_t e = gllel(eg) - Pr_t=0.91 - mu_t=rans_komg_mut(ix,iy,iz,e) + Pr_t=rans_turbPrandtl + mu_t=rans_mut(ix,iy,iz,e) - utrans = cpfld(ifield,2) if(ifield.eq.1) then + t(ix,iy,iz,e,4)=mu_t/cpfld(ifield,1) !store eddy viscosity for post processing udiff = cpfld(ifield,1)+mu_t + utrans = cpfld(ifield,2) elseif(ifield.eq.2) then udiff = cpfld(ifield,1)+mu_t*cpfld(ifield,2)/(Pr_t*cpfld(1,2)) - elseif(ifield.eq.3) then - udiff = cpfld(ifield,1)+rans_komg_mutsk(ix,iy,iz,e) - elseif(ifield.eq.4) then - udiff = cpfld(ifield,1)+rans_komg_mutso(ix,iy,iz,e) + utrans = cpfld(ifield,2) + elseif(ifield.eq.3) then !use rho and mu from field 1 + udiff = cpfld(1,1)+rans_mutsk(ix,iy,iz,e) + utrans = cpfld(1,2) + elseif(ifield.eq.4) then !use rho and mu from field 1 + udiff = cpfld(1,1)+rans_mutso(ix,iy,iz,e) + utrans = cpfld(1,2) endif return @@ -72,14 +75,17 @@ c----------------------------------------------------------------------- include 'NEKUSE' integer ix,iy,iz,e,eg - real rans_komg_kSrc,rans_komg_omgSrc + real rans_kSrc,rans_omgSrc + real rans_kDiag,rans_omgDiag e = gllel(eg) if(ifield.eq.3) then - qvol = rans_komg_kSrc(ix,iy,iz,e) + qvol = rans_kSrc(ix,iy,iz,e) + avol = rans_kDiag(ix,iy,iz,e) elseif(ifield.eq.4) then - qvol = rans_komg_omgSrc(ix,iy,iz,e) + qvol = rans_omgSrc(ix,iy,iz,e) + avol = rans_omgDiag(ix,iy,iz,e) else qvol = 0.0 endif @@ -112,30 +118,17 @@ c----------------------------------------------------------------------- include 'TOTAL' include 'NEKUSE' - real wd - common /walldist/ wd(lx1,ly1,lz1,lelv) - integer ix,iy,iz,e,eg - real darcy,utau,sigm,kmax,omax,yplus,fact,Re e = gllel(eg) - Re = 1.0/param(2) - darcy = 0.316/(Re**0.25) - utau = sqrt(darcy/8.0) - sigm = 0.6 - kmax = 4.5*utau*utau - omax = 0.5*utau*utau*Re - yplus = max(wd(ix,iy,iz,e)*utau*Re,1.0e-3) - - ux = 3.0/2.0*(1.0-y*y) + ux = 1.0 uy = 0.0 uz = 0.0 temp = 0.0 - fact = exp(-(log10(yplus)-1.0)**2/(2.0*sigm**2)) - if(ifield.eq.3) temp = kmax*fact - if(ifield.eq.4) temp = omax*fact + if(ifield.eq.3) temp = 0.01 + if(ifield.eq.4) temp = 0.2 return end @@ -145,23 +138,20 @@ c----------------------------------------------------------------------- include 'SIZE' include 'TOTAL' - real wd - common /walldist/ wd(lx1,ly1,lz1,lelv) - - real ypmin,ypmax,ypave,utmin,utmax,utave - - if (mod(istep,10).eq.0) then - call print_limits !monitor the solution - call y_p_limits(wd,ypmin,ypmax,ypave,utmin,utmax,utave) - if(nio.eq.0) then - write(6,256)'y_p+',ypmin,ypmax,ypave - write(6,256)'u_tau',utmin,utmax,utave - write(6,*) - endif + real xmax,xmin,ymax,ymin,zmax,zmin + real ptA(3),ptB(3) + + if(istep.eq.nsteps) then + call domain_size(xmin,xmax,ymin,ymax,zmin,zmax) + ptA(1) = 0.5*(xmin+xmax) + ptA(2) = ymax + ptA(3) = 0.5*(zmin+zmax) + ptB(1) = ptA(1) + ptB(2) = ymin + ptB(3) = ptA(3) + call RANSplot(ptA,ptB,1001) endif - 256 format(a15,3es13.4) - return end c----------------------------------------------------------------------- @@ -181,6 +171,9 @@ C enforce constant average velocity param(54) = -1 param(55) = 1.0 +c suppress runtime statistics + param(120) = nsteps + return end c----------------------------------------------------------------------- @@ -192,12 +185,18 @@ c----------------------------------------------------------------------- real wd common /walldist/ wd(lx1,ly1,lz1,lelv) - integer i,n + integer n,iel,ifc,id_face real xmin,xmax,ymin,ymax,scaley,scalex real glmin,glmax + integer ifld_k,ifld_t,m_id,w_id + real coeffs(30) !array for passing custom coeffs to RANS model + logical ifcoeffs !flag to use custom or default coeffs + n=nx1*ny1*nz1*nelv + +C rescale the domain BEFORE calling rans_init xmin=glmin(xm1,n) xmax=glmax(xm1,n) ymin=glmin(ym1,n) @@ -209,336 +208,208 @@ c----------------------------------------------------------------------- call cmult(xm1,scalex,n) call cmult(ym1,scaley,n) -c calculate wall distance for regularized k-omega - call rone(wd,n) - call sub2(wd,ym1,n) - return - end -c----------------------------------------------------------------------- - subroutine usrdat3() - implicit none - include 'SIZE' - include 'TOTAL' - - real wd - common /walldist/ wd(lx1,ly1,lz1,lelv) - - integer ifld_k,ifld_omg,m_id,w_id - real coeffs(21) !array for passing your own coeffs - logical ifcoeffs - - ifld_k = 3 !address of tke equation in t array - ifld_omg = 4 !address of omega equation in t array - ifcoeffs=.false. !set to true to pass your own coeffs - -C Supported models: -c m_id = 0 !regularized high-Re k-omega (no wall functions) - m_id = 1 !regularized low-Re k-omega -c m_id = 2 !regularized high-Re k-omega SST (no wall functions) -c m_id = 3 !regularized low-Re k-omega SST +C set BCs for velocity and temperature for 3rd party mesh +C tke and omega/tau are handled by rans_init +C do this BEFORE calling rans_init + +c do iel=1,nelv +c do ifc=1,2*ndim +c id_face = bc(5,ifc,iel,1) +c if (id_face.eq.1) then ! dirichlet (inlet) BCs +c cbc(ifc,iel,1) = 'v ' +c cbc(ifc,iel,2) = 't ' +c elseif (id_face.eq.2) then ! Wall / heat flux BCs +c cbc(ifc,iel,1) = 'W ' +c cbc(ifc,iel,2) = 'f ' +c elseif (id_face.eq.3) then ! Outlet BCs +c cbc(ifc,iel,1) = 'O ' +c cbc(ifc,iel,2) = 'I ' +c endif +c enddo +c enddo + + +C Setup RANS model, this MUST be done in usrdat2 + + ifld_k = 3 !field number for tke, t(1,1,1,1,ifld_k-1) + ifld_t = 4 !field number for omega/tau, t(1,1,1,1,ifld_t-1) + ifcoeffs = .false. !set to true to pass custom coefficients + +C Available models: +c m_id = 0 !regularized standard k-omega +c m_id = 1 !regularized low-Re k-omega +c m_id = 2 !regularized standard k-omega SST +c m_id = 3 !non-regularized standard k-omega (NOT SUPPORTED) + m_id = 4 !standard k-tau +c m_id = 5 !low-Re k-tau C Wall distance function: - w_id = 0 ! user specified +c w_id = 0 ! user specified c w_id = 1 ! cheap_dist (path to wall, may work better for periodic boundaries) -c w_id = 2 ! distf (coordinate difference, provides smoother function) + w_id = 2 ! distf (coordinate difference, provides smoother function) - call rans_komg_init(ifld_k,ifld_omg,ifcoeffs,coeffs,w_id,wd,m_id) + call rans_init(ifld_k,ifld_t,ifcoeffs,coeffs,w_id,wd,m_id) return end -C----------------------------------------------------------------------- - subroutine get_limits(phi,phimin,phimax,phiave,dphi,phip,rmsphi,n) +c----------------------------------------------------------------------- + subroutine usrdat3() implicit none include 'SIZE' include 'TOTAL' - integer i,n,ntot,iglsum - real phi(1),phip(1),phimin,phimax,phiave,dphi,rmsphi - real glmin,glmax,glsc2,glsum - - ntot=iglsum(n,1) - - rmsphi=0.0 - dphi=0.0 - if(istep.ge.1) then - do i=1,n - dphi=max(dphi,abs(phip(i)-phi(i))) - rmsphi=rmsphi+(phip(i)-phi(i))**2 - enddo - rmsphi=glsum(rmsphi,1) - rmsphi=sqrt(rmsphi/DBLE(ntot))/dt - endif - dphi=glmax(dphi,1) - dphi=dphi/dt - - phimin=glmin(phi,n) - phimax=glmax(phi,n) - phiave=glsc2(phi,bm1,n)/volvm1 - return end C----------------------------------------------------------------------- - subroutine get_limits_nodt(phi,phimin,phimax,phiave) - implicit none - include 'SIZE' - include 'TOTAL' - - integer n - real phi(1),phimin,phimax,phiave - real glmin,glmax,glsc2 - - n=nx1*ny1*nz1*nelv - - phimin=glmin(phi,n) - phimax=glmax(phi,n) - phiave=glsc2(phi,bm1,n)/volvm1 - - return - end -c----------------------------------------------------------------------- - subroutine y_p_limits(wd,ypmin,ypmax,ypave,utmin,utmax,utave) + subroutine RANSplot(pt1,pt2,lpts) implicit none include 'SIZE' include 'TOTAL' -C calculate min, max, and average y_p+ and u_tau values - - integer e,i,i0,i1,j,j0,j1,k,k0,k1,iw,jw,kw,i2,j2 - integer ipt,wpt,estrd,isd,jsd - real msk(lx1,ly1,lz1,lelv) - real gradu(lx1*ly1*lz1,3,3),wd(1) - real tau(3),norm(3),vsca,tauw,utau,rho,mu - real ypmin,ypmax,yp,ypave,vol,utmin,utmax,utave - real glmin,glmax,glsum - logical ifgrad, ifdid - - data ifdid /.false./ - save ifdid, msk - - ypmin=1.0d30 - ypmax=-1.0d30 - ypave=0.0 - utmin=1.0d30 - utmax=-1.0d30 - utave=0.0 - vol=0.0 - -C first build the mask - if(.not.ifdid)then - ifdid=.true. - call rone(msk,nx1*ny1*nz1*nelv) - do e=1,nelv - do isd=1,2*ndim - if(cbc(isd,e,1).eq.'W ') then - call backpts(i0,i1,j0,j1,k0,k1,isd) - do k=k0,k1 - do j=j0,j1 - do i=i0,i1 - msk(i,j,k,e)=0.0 - enddo - enddo - enddo - endif - enddo - do isd=1,2*ndim - if(cbc(isd,e,1).eq.'W ') then - call facind(i0,i1,j0,j1,k0,k1,lx1,ly1,lz1,isd) - do k=k0,k1 - do j=j0,j1 - do i=i0,i1 - msk(i,j,k,e)=1.0 - enddo - enddo - enddo - endif - enddo - enddo - call dssum(msk,nx1,ny1,nz1) !for elements with edges but not faces along a wall + real pt1(ldim),pt2(ldim) + integer npts,lpts,iplot + + character*32 fname + character*14 afmt + character*10 rfmt + integer intp_h,i,j,nt,nfld + save intp_h + logical ifset,ifdo + real dx,pts(lhis,ldim) + real fwrk(lx1*ly1*lz1*lelt,ldim+1+ldimt) + real fpts(lhis*(ldim+1+ldimt)) + real uout(lhis),vout(lhis),wout(lhis) + real prout(lhis),tout(lhis,ldimt) + character*4 outname(ldim+1+ldimt) + + real rwrk(lhis,ldim+1) + integer iwrk(lhis,3) + save rwrk,iwrk + + save ifdo,ifset + data ifdo /.true./ + data ifset /.true./ + + save iplot + data iplot /1/ + + if(.not.ifdo) return + + nt=lx1*ly1*lz1*nelt + + npts=max(lpts,2) + if(npts.gt.lhis) then + if(nio.eq.0) write(*,*) + & "Error, recompile with lhis in SIZE >= ",npts + ifdo=.false. + return endif - do e=1,nelv - ifgrad=.true. - do isd=1,2*ndim - if(cbc(isd,e,1).eq.'W ')then - estrd=(e-1)*nx1*ny1*nz1 - if(ifgrad)then - call gradm11(gradu(1,1,1),gradu(1,1,2),gradu(1,1,3),vx,e) - call gradm11(gradu(1,2,1),gradu(1,2,2),gradu(1,2,3),vy,e) - if(if3d) - & call gradm11(gradu(1,3,1),gradu(1,3,2),gradu(1,3,3),vz,e) - ifgrad=.false. - endif - call backpts(i0,i1,j0,j1,k0,k1,isd) - do k=k0,k1 - do j=j0,j1 - do i=i0,i1 - if(msk(i,j,k,e).lt.0.5) then - iw=i - jw=j - kw=k - if (isd.eq.1) then - jw=1 - elseif(isd.eq.2) then - iw=nx1 - elseif(isd.eq.3) then - jw=ny1 - elseif(isd.eq.4) then - iw=1 - elseif(isd.eq.5) then - kw=1 - else - kw=nx1 - endif - call getSnormal(norm,iw,jw,kw,isd,e) - ipt=i +(j -1)*nx1+(k -1)*nx1*ny1 - wpt=iw+(jw-1)*nx1+(kw-1)*nx1*ny1 - - mu=vdiff(iw,jw,kw,e,1) - rho=vtrans(iw,jw,kw,e,1) - - do i2=1,ldim - tau(i2)=0.0 - do j2=1,ldim - tau(i2)=tau(i2)+ - & mu*(gradu(wpt,i2,j2)+gradu(wpt,j2,i2))*norm(j2) - enddo - enddo - - vsca=0.0 - do i2=1,ldim - vsca=vsca+tau(i2)*norm(i2) - enddo - - tauw=0.0 - do i2=1,ldim - tauw=tauw+(tau(i2)-vsca*norm(i2))**2 - enddo - tauw=sqrt(tauw) - utau=sqrt(tauw/rho) - yp=wd(ipt+estrd)*utau*rho/mu - ypmin=min(ypmin,yp) - ypmax=max(ypmax,yp) - ypave=ypave+yp*bm1(i,j,k,e) - utmin=min(utau,utmin) - utmax=max(utau,utmax) - utave=utave+utau*bm1(i,j,k,e) - vol=vol+bm1(i,j,k,e) - endif - enddo - enddo - enddo - endif + call rzero(pts,npts*ndim) + do j=1,ndim + pts(1,j)=pt1(j) + dx=(pt2(j)-pt1(j))/(real(npts-1)) + do i=2,npts + pts(i,j)=pts(i-1,j)+dx enddo enddo - ypmin=glmin(ypmin,1) - ypmax=glmax(ypmax,1) - ypave=glsum(ypave,1) - utmin=glmin(utmin,1) - utmax=glmax(utmax,1) - utave=glsum(utave,1) - vol=glsum(vol,1) - ypave=ypave/vol - utave=utave/vol - - return - end -c----------------------------------------------------------------------- - subroutine backpts(i0,i1,j0,j1,k0,k1,isd) - implicit none - include 'SIZE' - - integer i0,i1,j0,j1,k0,k1,isd - - i0=1 - j0=1 - k0=1 - i1=nx1 - j1=ny1 - k1=nz1 - if(isd.eq.1) then - j0=2 - j1=2 - elseif(isd.eq.2) then - i0=nx1-1 - i1=nx1-1 - elseif(isd.eq.3) then - j0=ny1-1 - j1=ny1-1 - elseif(isd.eq.4) then - i0=2 - i1=2 - elseif(isd.eq.5) then - k0=2 - k1=2 - elseif(isd.eq.6) then - k0=nz1-1 - k1=nz1-1 + if(ifset)then + ifset=.false. + call interp_setup(intp_h,0.0,0,nelt) endif - return - end -c----------------------------------------------------------------------- - subroutine print_limits - implicit none - include 'SIZE' - include 'TOTAL' - - real vol,glsum,glmin,glmax,glsc2 - real tmp(lx1*ly1*lz1*lelv) - real prlagl(lx2*ly2*lz2*lelv) - integer i,n1,n2,nt - - real uxmin,uxmax,uxave,uymin,uymax,uyave,uzmin,uzmax,uzave - real prmin,prmax,prave - real thmin(ldimt),thmax(ldimt),thave(ldimt) - real rmsux,rmsuy,rmsuz,rmspr,rmsth(ldimt) - real dux,duy,duz,dpr,dth(ldimt) - character*15 tname - - save prlagl - - n1=nx1*ny1*nz1*nelv - n2=nx2*ny2*nz2*nelv - nt=nx1*ny1*nz1*nelt - - call get_limits(vx,uxmin,uxmax,uxave,dux,vxlag,rmsux,n1) - call get_limits(vy,uymin,uymax,uyave,duy,vylag,rmsuy,n1) - if(if3d) call get_limits(vz,uzmin,uzmax,uzave,duz,vzlag,rmsuz,n1) - call get_limits(pr,prmin,prmax,prave,dpr,prlagl,rmspr,n2) - if(ifheat)then - do i=1,npscal+1 - if(idpss(i).eq.0) call get_limits(t(1,1,1,1,i),thmin(i) !Helmholtz solver - & ,thmax(i),thave(i),dth(i),tlag(1,1,1,1,1,i),rmsth(i),nt) - if(idpss(i).eq.1) call get_limits_nodt(t(1,1,1,1,i) !CVODE solver - & ,thmin(i),thmax(i),thave(i)) + nfld=0 + if(ifvo) then + write(outname(1),'(a4)')"VELX" + write(outname(2),'(a4)')"VELY" + call copy(fwrk(1,1),vx,nt) + call copy(fwrk(1,2),vy,nt) + nfld=2 + endif + if(if3d.and.ifvo)then + nfld=nfld+1 + write(outname(nfld),'(a4)')"VELZ" + call copy(fwrk(1,nfld),vz,nt) + endif + if(ifpo) then + nfld=nfld+1 + write(outname(nfld),'(a4)')"PRES" + call copy(fwrk(1,nfld),pr,nt) + endif + if(ifheat) then + if(ifto) then + nfld=nfld+1 + write(outname(nfld),'(a4)')"TEMP" + call copy(fwrk(1,nfld),t,nt) + endif + do i=1,ldimt-1 + if(ifpsco(i)) then + nfld=nfld+1 + write(outname(nfld),'(a2,i2)')"PS",i + call copy(fwrk(1,nfld),t(1,1,1,1,i+1),nt) + endif enddo endif - if(nio.eq.0) then - write(*,*) - write(*,254) 'limits','min','max','ave','max d/dt','rms d/dt' - write(*,255) 'u velocity',uxmin,uxmax,uxave,dux,rmsux - write(*,255) 'v velocity',uymin,uymax,uyave,duy,rmsuy - if(if3d) write(*,255) 'w velocity',uzmin,uzmax,uzave,duz,rmsuz - write(*,255) 'pressure',prmin,prmax,prave,dpr,rmspr - if(ifheat) then - do i=1,npscal+1 - if(i.eq.1) write(tname,'(a15)') "temperature" - if(i.gt.1) write(tname,'(a14,i1)') "PS ",i-1 - if(idpss(i).eq.0)write(*,255) - & tname,thmin(i),thmax(i),thave(i),dth(i),rmsth(i) - if(idpss(i).eq.1)write(*,256) - & tname,thmin(i),thmax(i),thave(i),'--','--' + if(nfld.gt.0) then + call blank(fname,32) + if(iplot.lt.10) then + write(fname,'(a,i1,a)') "plot",iplot,".dat" + elseif(iplot.lt.100) then + write(fname,'(a,i2,a)') "plot",iplot,".dat" + else + write(fname,'(a,i3,a)') "plot",iplot,".dat" + endif + + if(nio.eq.0) then + write(*,*)' Writing line plot data to file ',fname + if(if3d)then + write(*,'(7x,3es15.6)')pt1(1),pt1(2),pt1(3) + write(*,'(7x,3es15.6)')pt2(1),pt2(2),pt2(3) + else + write(*,'(7x,2es15.6)')pt1(1),pt1(2) + write(*,'(7x,2es15.6)')pt2(1),pt2(2) + endif + write(*,*) + endif + + call interp_nfld(fpts,fwrk,nfld,pts(1,1),pts(1,2),pts(1,3),npts + & ,iwrk,rwrk,lhis,.true.,intp_h) + + call blank(afmt,14) + call blank(rfmt,10) + if(if3d) then + write(afmt,'(a1,i2,a11)')"(",nfld+3,"a16,es16.8)" + write(rfmt,'(a1,i2,a7)')"(",nfld+3,"es16.8)" + else + write(afmt,'(a1,i2,a11)')"(",nfld+2,"a16,es16.8)" + write(rfmt,'(a1,i2,a7)')"(",nfld+2,"es16.8)" + endif + + if(nio.eq.0) then + open(unit=10,file=fname,status='unknown',form='formatted') + if(if3d) then + write(10,afmt)"X","Y","Z",(outname(i),i=1,nfld),time + else + write(10,afmt)"X","Y",(outname(i),i=1,nfld),time + endif + do i=1,npts + if(if3d) then + write(10,rfmt)pts(i,1),pts(i,2),pts(i,3) + & ,(fpts(i+j),j=0,(npts*nfld-1),npts) + else + write(10,rfmt)pts(i,1),pts(i,2) + & ,(fpts(i+j),j=0,(npts*nfld-1),npts) + endif enddo endif - write(*,*) - endif - 254 format(a15,5a13) - 255 format(a15,5es13.4) - 256 format(a15,3es13.4,2a13) + close(10) + + iplot=iplot+1 + endif return end diff --git a/RANSChannel/README.md b/RANSChannel/README.md index f63168f..95a8bd7 100644 --- a/RANSChannel/README.md +++ b/RANSChannel/README.md @@ -1,3 +1,3 @@ # RANS modle for 2D chnnel flow -This example demonstrates use of kOmg turbulence model. +This example demonstrates the use of the standar k-tau turbulence model. diff --git a/RANSChannel/SIZE b/RANSChannel/SIZE index d2b48a8..3762a0c 100644 --- a/RANSChannel/SIZE +++ b/RANSChannel/SIZE @@ -14,25 +14,26 @@ c parameter (lxd=12) ! GL points for over-integration (dealiasing) parameter (lx2=lx1-0) ! GLL points for pressure (lx1 or lx1-2) - parameter (lelg=64) ! max number of global elements + parameter (lelg=30) ! max number of global elements parameter (lpmin=1) ! min number of MPI ranks parameter (lelt=lelg/lpmin + 3) ! max number of local elements per MPI rank - parameter (ldimt=3) ! max auxiliary fields (temperature + scalars) + parameter (ldimt=4) ! max auxiliary fields (temperature + scalars) ! OPTIONAL parameter (ldimt_proj=1) ! max auxiliary fields residual projection parameter (lelr=lelt) ! max number of local elements per restart file - parameter (lhis=1000) ! max history/monitoring points + parameter (lhis=1001) ! max history/monitoring points parameter (maxobj=1) ! max number of objects parameter (lpert=1) ! max number of perturbations parameter (toteq=1) ! max number of conserved scalars in CMT parameter (nsessmax=1) ! max sessions to NEKNEK parameter (lxo=lx1) ! max GLL points on output (lxo>=lx1) - parameter (mxprev=20,lgmres=30) ! max dim of projection & Krylov space + parameter (mxprev=20) ! max dim of projection space + parameter (lgmres=30) ! max dim Krylov space parameter (lorder=3) ! max order in time - parameter (lx1m=lx1) ! GLL points mesh solver + parameter (lx1m=lx1) ! GLL points mesh solver parameter (lfdm=0) ! unused - parameter (lelx=3,lely=8,lelz=1) ! global tensor mesh dimensions + parameter (lelx=1,lely=1,lelz=1) ! global tensor mesh dimensions parameter (lbelt=1) ! lelt for mhd parameter (lpelt=1) ! lelt for linear stability From 4da79dd8b61dcbe78dc558c474cb1557a99ffbb0 Mon Sep 17 00:00:00 2001 From: Dillon Shaver Date: Tue, 6 Apr 2021 09:34:59 -0500 Subject: [PATCH 3/5] Instead of bc(5,ifc,iel,1) use boundaryID(ifc,iel) --- RANSChannel/RANSchan.usr | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RANSChannel/RANSchan.usr b/RANSChannel/RANSchan.usr index d1e694e..a10103a 100644 --- a/RANSChannel/RANSchan.usr +++ b/RANSChannel/RANSchan.usr @@ -215,7 +215,7 @@ C do this BEFORE calling rans_init c do iel=1,nelv c do ifc=1,2*ndim -c id_face = bc(5,ifc,iel,1) +c id_face = boundaryID(ifc,iel) c if (id_face.eq.1) then ! dirichlet (inlet) BCs c cbc(ifc,iel,1) = 'v ' c cbc(ifc,iel,2) = 't ' From b0e5199e0e9bab6d4cb6f7105487e8dc00242f16 Mon Sep 17 00:00:00 2001 From: Dillon Shaver Date: Mon, 24 May 2021 16:07:26 -0500 Subject: [PATCH 4/5] fix typos in RANSchannel --- RANSChannel/RANSchan.usr | 2 +- RANSChannel/README.md | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/RANSChannel/RANSchan.usr b/RANSChannel/RANSchan.usr index a10103a..3bf734e 100644 --- a/RANSChannel/RANSchan.usr +++ b/RANSChannel/RANSchan.usr @@ -27,7 +27,7 @@ c----------------------------------------------------------------------- e = gllel(eg) - Pr_t=rans_turbPrandtl + Pr_t=rans_turbPrandtl() mu_t=rans_mut(ix,iy,iz,e) if(ifield.eq.1) then diff --git a/RANSChannel/README.md b/RANSChannel/README.md index 95a8bd7..2ef2404 100644 --- a/RANSChannel/README.md +++ b/RANSChannel/README.md @@ -1,3 +1,3 @@ -# RANS modle for 2D chnnel flow +# RANS model for 2D channel flow -This example demonstrates the use of the standar k-tau turbulence model. +This example demonstrates the use of the standard k-tau turbulence model. From 8da083a091a86f5bd3847c821721ebfa3382b471 Mon Sep 17 00:00:00 2001 From: Dillon Shaver Date: Wed, 6 May 2026 13:34:25 -0500 Subject: [PATCH 5/5] fix typo in prepost call --- turbChannel/turbChannel.usr | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/turbChannel/turbChannel.usr b/turbChannel/turbChannel.usr index d0647a2..653101c 100644 --- a/turbChannel/turbChannel.usr +++ b/turbChannel/turbChannel.usr @@ -119,7 +119,7 @@ c----------------------------------------------------------------------- call create_obj(iobj_wall,bIDs,1) nm = iglsum(nmember(iobj_wall),1) if(nid.eq.0) write(6,*) 'obj_wall nmem:', nm - call prepost(.true.,' ') + call prepost(.true.,' ') endif ubar = glsc2(vx,bm1,n)/volvm1