Permalink
Browse files

Some corrections to the calculation of the dissipation rates by

calculating derivatives using  neighboring grid points. This version
now gives converged results for the kinetic and thermal dissipation
rates also on very low grid resolutions. Tested on 64^3, 96^3, 128^3,
and 192^3. Observed  convergence well within 1% when averaged over 300
dimensionless time units.
  • Loading branch information...
stevensrjam committed Feb 17, 2016
1 parent b54ebe6 commit 2fc64fbe8bdf46d6bb8b83c77e7777f52410a520
Showing with 51 additions and 63 deletions.
  1. +51 −63 CalcDissipationNu.F90
View
@@ -18,115 +18,103 @@ subroutine CalcDissipationNu
implicit none
integer :: i,j,k
integer :: imm,ipp,jmm,jpp,kmm,kpp,kp
real :: udx3_m,udx3_c
integer :: ip,jp,kp
real :: dxm(nxm),dxc(nxm),lem(nxm),lec(nxm)
real :: hxx,hxy,hxz,hyx,hyy,hyz,hzx,hzy,hzz
real :: tx,ty,tz
real :: nu_th,nu_mu,volt
real :: udz,udy,dissipte,dissipth
real :: dissipte,dissipte2,dissipth
nu_mu = 0.0d0
nu_th = 0.0d0
udy=dy
udz=dz
! RS: Calculated geometrical quantities outside the main loop.
! RS: Could potentially be moved to initialization routine
do k=1,nxm
kp=k+1
dxm(k)=1.0/(xm(kp)-xm(k))
dxc(k)=1.0/(xc(kp)-xc(k))
lem(k)=(xc(kp)-xc(k))*dx
lec(k)=(xm(kp)-xm(k))*dx
enddo
lec(nxm)=(xm(nxm)-xm(nxm-1))*dx
dxm(nxm)=1.0/(xm(nxm)-xm(nxm-1))
!$OMP PARALLEL DO &
!$OMP DEFAULT(none) &
!$OMP SHARED(xstart,xend,g3rm,pra,vz,vy,vx,temp) &
!$OMP SHARED(udz,udy,udx3m,udx3c) &
!$OMP SHARED(nzm,nym,nxm,ren,pec) &
!$OMP SHARED(kpv,kmv,xc,xm) &
!$OMP SHARED(xstart,xend,vz,vy,vx,temp) &
!$OMP SHARED(nxm,ren,pec,pra) &
!$OMP SHARED(dxc,dxm,lec,lem) &
!$OMP SHARED(disste,dissth) &
!$OMP PRIVATE(i,j,k,imm,ipp,jmm,jpp,kp,kpp,kmm) &
!$OMP PRIVATE(udx3_m,udx3_c,dissipte,dissipth) &
!$OMP PRIVATE(i,j,k,ip,jp,kp) &
!$OMP PRIVATE(dissipte,dissipte2,dissipth) &
!$OMP PRIVATE(hxx,hxy,hxz,hyx,hyy,hyz,hzx,hzy,hzz) &
!$OMP PRIVATE(tx,ty,tz) &
!$OMP REDUCTION(+:nu_mu) &
!$OMP REDUCTION(+:nu_th)
do i=xstart(3),xend(3)
imm= i-1
ipp= i+1
ip= i+1
do j=xstart(2),xend(2)
jmm=j-1
jpp=j+1
jp=j+1
do k=1,nxm
kp=k+1
kpp=kpv(k)
kmm=kmv(k)
udx3_m=1.0/(xm(kpp)-xm(kmm))
udx3_c=1.0/(xc(kp)-xc(kmm))
!
! Viscous dissipation rate
!
! Viscous dissipation rate
! 1 | | 2
! ---- | nabla u|
! Re | |
!
hxx=(vx(kp,j,i)-vx(k,j,i))*udx3c(k)
hxy=(vx(k,jpp,i)-vx(k,j,i))*udy
hxz=(vx(k,j,ipp)-vx(k,j,i))*udz
hxx=(vx(kp,j,i)-vx(k,j,i))*dxc(k)
hxy=(vx(k,jp,i)-vx(k,j,i))*dy
hxz=(vx(k,j,ip)-vx(k,j,i))*dz
hyx=(vy(kpp,j,i)-vy(k,j,i))*udx3m(k)
hyy=(vy(k,jpp,i)-vy(k,j,i))*udy
hyz=(vy(k,j,ipp)-vy(k,j,i))*udz
hyx=(vy(kp,j,i)-vy(k,j,i))*dxm(k)
hyy=(vy(k,jp,i)-vy(k,j,i))*dy
hyz=(vy(k,j,ip)-vy(k,j,i))*dz
hzx=(vz(kpp,j,i)-vz(k,j,i))*udx3m(k)
hzy=(vz(k,jpp,i)-vz(k,j,i))*udy
hzz=(vz(k,j,ipp)-vz(k,j,i))*udz
hzx=(vz(kp,j,i)-vz(k,j,i))*dxm(k)
hzy=(vz(k,jp,i)-vz(k,j,i))*dy
hzz=(vz(k,j,ip)-vz(k,j,i))*dz
dissipte = 2.0*(hxx**2+hyy**2+hzz**2)+ &
(hyz+hzy)**2+(hxz+hzx)**2+(hxy+hyx)**2
dissipte = (hxx*hxx+hxy*hxy+hxz*hxz)
dissipte2 = (hyx*hyx+hyy*hyy+hyz*hyz)+(hzx*hzx+hzy*hzy+hzz*hzz)
nu_mu = nu_mu + dissipte*lem(k)*pra+dissipte2*lec(k)*pra
nu_mu = nu_mu + dissipte*g3rm(k)*pra
!
! Thermal gradient dissipation rate
!
! 1 | | 2
! ---- | nabla T|
! Pe | |
!
tx=(temp(kp,j,i )-temp(kmm,j,i))*udx3_c
ty=(temp(k,jpp,i)-temp(k,jmm,i))*udy*0.5
tz=(temp(k,j,ipp)-temp(k,j,imm))*udz*0.5
tx=(temp(kp,j,i)-temp(k,j,i))*dxc(k)
ty=(temp(k,jp,i)-temp(k,j,i))*dy
tz=(temp(k,j,ip)-temp(k,j,i))*dz
dissipth = tx*tx + ty*ty + tz*tz
nu_th = nu_th+dissipth*g3rm(k)
nu_th = nu_th+dissipth*lem(k)
!$OMP CRITICAL
disste(k) = disste(k) + dissipte / (ren*real(nym)*real(nzm))
dissth(k) = dissth(k) + dissipth / (pec*real(nym)*real(nzm))
disste(k) = disste(k) + (dissipte + dissipte2) / (ren*real(nym)*real(nzm))
dissth(k) = dissth(k) + dissipth / (pec*real(nym)*real(nzm))
!$OMP END CRITICAL
end do
end do
end do
!$OMP END PARALLEL DO
call MpiSumRealScalar(nu_th)
call MpiSumRealScalar(nu_mu)
volt = 1.d0/(real(nxm)*real(nzm)*real(nym))
if(ismaster) then
nu_mu = nu_mu*volt + 1
nu_th = nu_th*volt
open(92,file='nu_diss.out',status='unknown',access='sequential', &
position='append')
write(92,*) time,nu_mu,nu_th
close(92)
endif
return
end
if(ismaster) then
nu_mu = nu_mu*volt + 1
nu_th = nu_th*volt
open(92,file='nu_diss.out',status='unknown',access='sequential',position='append')
write(92,*) time,nu_mu,nu_th
close(92)
endif
return
end

0 comments on commit 2fc64fb

Please sign in to comment.