Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • GOzler-master-patch-22703
  • master
  • v1.0.0
  • v1.1.0
4 results

Target

Select target project
  • Holger/pafix
1 result
Select Git revision
  • GOzler-master-patch-22703
  • master
  • v1.0.0
  • v1.1.0
4 results
Show changes
Commits on Source (8)
Hide whitespace changes
Inline Side-by-side
Showing
with 228 additions and 113 deletions
clean.sh 100644 → 100755
View file @ f40b4e22
File mode changed from 100644 to 100755
documentation/doxygen_codeGuide/codeGuide 100644 → 100755
View file @ f40b4e22
File mode changed from 100644 to 100755
documentation/doxygen_codeGuide/pafiXlogo.png 100644 → 100755
View file @ f40b4e22
File mode changed from 100644 to 100755
input/input_example.dat 100644 → 100755
View file @ f40b4e22
......@@ -12,6 +12,8 @@ start time-step (-), time-steps to compute (-) =
1,10
interval of writing result files, restart files (-,-) =
10,10000
turbulence model [-/Smagorinsky] =
Smagorinsky
CFL number (-), max number interations (-), relative tolerance (-) =
0.4,200,5.E-3
in/initial flow (m/s), pressure gradient (N/m**3) =
......
monitor.plt 100644 → 100755
View file @ f40b4e22
File mode changed from 100644 to 100755
run.sh 100644 → 100755
View file @ f40b4e22
File mode changed from 100644 to 100755
src/COPYING 100644 → 100755
View file @ f40b4e22
File mode changed from 100644 to 100755
src/bc.f90 100644 → 100755
View file @ f40b4e22
......@@ -10,6 +10,7 @@
wval=0._pr
call sync(myu); call sync(myv); call sync(myw) ! between procs through the sync routines
if (bcx.eq.'w') then ! x:
if (myid.eq.0) then
do 1 m=1,gc
......@@ -26,9 +27,9 @@
2 enddo
endif
elseif (bcx.eq.'p') then
call sync(u); call sync(v); call sync(w) ! done through the sync routines
! between procs through the sync routines
elseif (bcx.eq.'i') then
! see inflow routine
! see inflow routine
if (myid.eq.nrprocs-1) then
do 8 m=1,gc
myu(imax+m-1,:,:)=myu(imax-1,:,:)
......@@ -85,7 +86,6 @@
!> @brief inflow boundary condition
subroutine inflow
use var
use parallel
real(kind=pr),dimension(dimj*diml) :: random
real(kind=pr) :: alpha,dist
integer :: i,j,l,m,mm
......@@ -116,6 +116,7 @@
real(kind=pr) :: myvar(ii,jj,ll)
integer :: m
call sync(myvar) ! between procs through the sync routines
if (bcx.eq.'w') then
if (myid.eq.0) then
do 1 m=1,gc
......@@ -128,7 +129,7 @@
2 enddo
endif
elseif (bcx.eq.'p') then
call sync(myvar) ! done through the sync routines
! between procs through the sync routines
elseif (bcx.eq.'i') then
if (myid.eq.0) then
do 7 m=1,gc
......@@ -178,6 +179,7 @@
real(kind=pr) :: wval
integer :: m
call sync(myvar) ! between procs through the sync routines
if (bcx.eq.'w') then
if (myid.eq.0) then
do 1 m=1,gc
......@@ -190,7 +192,7 @@
2 enddo
endif
elseif (bcx.eq.'p') then
call sync(myvar) ! done through the sync routines
! between procs through the sync routines
elseif (bcx.eq.'i') then
if (myid.eq.0) then
do 7 m=1,gc
......@@ -243,9 +245,7 @@
real(kind=pr) :: myvar(ii,jj,ll)
integer :: m
if (bcx.eq.'p') then
call syncLEsource(myvar) ! done through the sync routines
endif
call syncLEsource(myvar) ! between procs through the sync routines
if (bcy.eq.'p') then
do 4 m=1,gc
......
src/electrostatics.f90 100644 → 100755
View file @ f40b4e22
......@@ -101,7 +101,7 @@
!####################################################################
!> @author Holger Grosshans
!> @brief compute electric forces on particles from E-field
!> @brief compute specific electric forces on particles from E-field
!> at seeding time-step E is not available, only Coulomb forces
subroutine forcesGauss
use var
......@@ -126,7 +126,7 @@
!####################################################################
!> @author Holger Grosshans
!> @brief compute electric forces on particles from Coulomb's law
!> @brief compute specific electric forces on particles from Coulomb's law
subroutine forcesCoulomb
use var
use parallel
......
src/fluid.f90 100644 → 100755
View file @ f40b4e22
......@@ -5,14 +5,14 @@
subroutine solveFluid
use var
use parallel
real(kind=pr),dimension(ii,jj,ll) :: du,dv,dw,dp,ddp,massRes
real(kind=pr),dimension(ii,jj,ll) :: du,dv,dw,dp,massRes
real(kind=pr) :: err(4),err0(4),errL2
integer itout,it,itmommax,itpCorr1max,itpCorr2max,velCorr1,velCorr2
!> itmax: outer iterations
itmommax=1; itpCorr1max=1; itpCorr2max=0 ! sequential SIMPLE
!itmommax=4; itpCorr1max=4; itpCorr2max=0 ! SIMPLE
!itmommax=4; itpCorr1max=4; itpCorr2max=4 ! PISO
!itmommax>1; itpCorr1max>1; itpCorr2max=0 ! SIMPLE
!itmommax>1; itpCorr1max>1; itpCorr2max>1 ! PISO
velCorr1=0; velCorr2=0 ! velocity correction yes/no (1/0)
if ((bcx.eq.'i').and.(myid.eq.0)) call inflow
......@@ -20,11 +20,12 @@
do 1 itout=1,itmax
if ((mod(it,10).eq.0).and.(npTot.ne.0)) call momentumCoupling
if (turbModel.eq.'Smagorinsky') call turbSmagorinsky
do 2 it=1,itmommax ! solve momentum eq
call momx(du); call momy(dv); call momz(dw)
call bcUVW(du,dv,dw)
u=u+du*urfu; v=v+dv*urfv; w=w+dw*urfw
u=u+du*urfu; v=v+dv*urfv; w=w+dw*urfw
call bcUVW(u,v,w)
2 enddo
call mass(massRes,u,v,w) ! first pressure correction
......@@ -35,10 +36,10 @@
3 enddo
p= p+dp*urfp
if (velCorr2.eq.1) then ! first velocity correction
if (velCorr1.eq.1) then ! first velocity correction
call velocityCorrection(du,dv,dw,dp)
call bcUVW(du,dv,dw)
u=u+du*urfu; v=v+dv*urfv; w=w+dw*urfw
u=u+du*urfu; v=v+dv*urfv; w=w+dw*urfw
endif
if (itpCorr2max.ge.1) then
......@@ -48,12 +49,12 @@
call pressureCorrection(dp,massRes)
call bcNeumann(dp)
4 enddo
p= p+dp*urfp
p=p+dp*urfp
if (velCorr2.eq.1) then ! second velocity correction
call velocityCorrection(du,dv,dw,dp)
call bcUVW(du,dv,dw)
u=u+du*urfu; v=v+dv*urfv; w=w+dw*urfw
u=u+du*urfu; v=v+dv*urfv; w=w+dw*urfw
endif
endif
......@@ -67,7 +68,7 @@
1 enddo
if (myid.eq.0) write(*,'(a,i3,3(a,es8.2e2))') 'fluid |L2 #',itout-1,' = ',errL2, &
' |mom. = ',sum(err(1:3)/err0(1:3))/3._pr,' |pc. = ',err(4)/err0(4)
' |mom. = ',sum(err(1:3)/err0(1:3))/3._pr,' |pc. = ',err(4)/err0(4)
end
......@@ -125,7 +126,7 @@
ux= (myu(i,j,l)-myu(i-1,j,l))/(xf(i)-xf(i-1)) ! 2nd order central
vy= (myv(i,j,l)-myv(i,j-1,l))/(yf(j)-yf(j-1))
wz= (myw(i,j,l)-myw(i,j,l-1))/(zf(l)-zf(l-1))
massRes(i,j,l)=-(ux+vy+wz)
massRes(i,j,l)= -(ux+vy+wz)
enddo; enddo; enddo
end
......@@ -136,8 +137,7 @@
subroutine momx(du)
use var
real(kind=pr),dimension(ii,jj,ll) :: tu,qu,du
real(kind=pr) :: ua,va,wa,px,uux,vuy,wuz, &
uxx,uyy,uzz,dudt,Cdudt,Cviscx,erru
real(kind=pr) :: ua,va,wa,px,uux,vuy,wuz,uxx,uyy,uzz,dudt,Cdudt,Cviscx
integer :: i,j,l
tu=0._pr; qu=0._pr; du=0._pr
......@@ -150,7 +150,7 @@
va= 0.25_pr*(v(i,j,l)+v(i,j-1,l)+v(i+1,j,l)+v(i+1,j-1,l))
wa= 0.25_pr*(w(i,j,l)+w(i,j,l-1)+w(i+1,j,l)+w(i+1,j,l-1))
! convective terns (2nd order central)
! convective terms (2nd order central)
uux= ua*(u(i+1,j,l)-u(i-1,j,l))/(xf(i+1)-xf(i-1))
vuy= va*(u(i,j+1,l)-u(i,j-1,l))/(yc(j+1)-yc(j-1))
wuz= wa*(u(i,j,l+1)-u(i,j,l-1))/(zc(l+1)-zc(l-1))
......@@ -171,13 +171,13 @@
(2._pr*tau01*dt)
! momentum eq -lhs+rhs (m/s**2)
tu(i,j,l)= - dudt -(uux+vuy+wuz) -(px+dpdx)/rhof + nuf*(uxx+uyy+uzz) + Fsx(i,j,l)
tu(i,j,l)= - dudt -(uux+vuy+wuz) -(px+dpdx)/rhof + (nuf+nuf_ru(i,j,l))*(uxx+uyy+uzz) + Fsx(i,j,l)
! coefficients of u(i,j,l) of momentum eq
Cdudt= (1._pr+tau01)/(2._pr*tau01*dt)
Cviscx= -2._pr*nuf * ( 1._pr/((xf(i+1)-xf(i))*(xf(i)-xf(i-1))) &
+ 1._pr/((yc(j+1)-yc(j))*(yc(j)-yc(j-1))) &
+ 1._pr/((zc(l+1)-zc(l))*(zc(l)-zc(l-1))) )
Cviscx= -2._pr*(nuf+nuf_ru(i,j,l)) * ( 1._pr/((xf(i+1)-xf(i))*(xf(i)-xf(i-1))) &
+ 1._pr/((yc(j+1)-yc(j))*(yc(j)-yc(j-1))) &
+ 1._pr/((zc(l+1)-zc(l))*(zc(l)-zc(l-1))) )
! coefficients lhs-rhs
qu(i,j,l)= Cdudt - Cviscx
......@@ -194,8 +194,7 @@
subroutine momy(dv)
use var
real(kind=pr),dimension(ii,jj,ll) :: tv,qv,dv
real(kind=pr) :: ua,va,wa,py,uvx,vvy,wvz, &
vxx,vyy,vzz,dvdt,Cdvdt,Cviscy
real(kind=pr) :: ua,va,wa,py,uvx,vvy,wvz,vxx,vyy,vzz,dvdt,Cdvdt,Cviscy
integer :: i,j,l
tv=0._pr; qv=0._pr; dv=0._pr
......@@ -208,7 +207,7 @@
va= v(i,j,l)
wa= 0.25_pr*(w(i,j,l)+w(i,j,l-1)+w(i,j+1,l)+w(i,j+1,l-1))
! convective terns (2nd order central)
! convective terms (2nd order central)
uvx= ua*(v(i+1,j,l)-v(i-1,j,l))/(xc(i+1)-xc(i-1))
vvy= va*(v(i,j+1,l)-v(i,j-1,l))/(yf(j+1)-yf(j-1))
wvz= wa*(v(i,j,l+1)-v(i,j,l-1))/(zc(l+1)-zc(l-1))
......@@ -229,13 +228,13 @@
(2._pr*tau01*dt)
! momentum eq. -lhs+rhs
tv(i,j,l)= - dvdt - (uvx+vvy+wvz) - py/rhof + nuf*(vxx+vyy+vzz) + Fsy(i,j,l)
tv(i,j,l)= - dvdt - (uvx+vvy+wvz) - py/rhof + (nuf+nuf_rv(i,j,l))*(vxx+vyy+vzz) + Fsy(i,j,l)
! coefficients of v(i,j,l) of momentum eq
Cdvdt= (1._pr+tau01)/(2._pr*tau01*dt)
Cviscy= -2._pr*nuf * ( 1._pr/((xc(i+1)-xc(i))*(xc(i)-xc(i-1))) &
+ 1._pr/((yf(j+1)-yf(j))*(yf(j)-yf(j-1))) &
+ 1._pr/((zc(l+1)-zc(l))*(zc(l)-zc(l-1))) )
Cviscy= -2._pr*(nuf+nuf_rv(i,j,l)) * ( 1._pr/((xc(i+1)-xc(i))*(xc(i)-xc(i-1))) &
+ 1._pr/((yf(j+1)-yf(j))*(yf(j)-yf(j-1))) &
+ 1._pr/((zc(l+1)-zc(l))*(zc(l)-zc(l-1))) )
! coefficients lhs-rhs
qv(i,j,l)= Cdvdt - Cviscy
......@@ -252,8 +251,7 @@
subroutine momz(dw)
use var
real(kind=pr),dimension(ii,jj,ll) :: tw,qw,dw
real(kind=pr) :: ua,va,wa,pz,uwx,vwy,wwz, &
wxx,wyy,wzz,dwdt,Cdwdt,Cviscz
real(kind=pr) :: ua,va,wa,pz,uwx,vwy,wwz,wxx,wyy,wzz,dwdt,Cdwdt,Cviscz
integer :: i,j,l
tw=0._pr; qw=0._pr; dw=0._pr
......@@ -266,7 +264,7 @@
va= 0.25_pr*(v(i,j,l)+v(i,j-1,l)+v(i,j,l+1)+v(i,j-1,l+1))
wa= w(i,j,l)
! convective terns (2nd order central)
! convective terms (2nd order central)
uwx= ua*(w(i+1,j,l)-w(i-1,j,l))/(xc(i+1)-xc(i-1))
vwy= va*(w(i,j+1,l)-w(i,j-1,l))/(yc(j+1)-yc(j-1))
wwz= wa*(w(i,j,l+1)-w(i,j,l-1))/(zf(l+1)-zf(l-1))
......@@ -287,13 +285,13 @@
(2._pr*tau01*dt)
! momentum eq. -lhs+rhs
tw(i,j,l)= - dwdt - (uwx+vwy+wwz) - pz/rhof + nuf*(wxx+wyy+wzz) + Fsz(i,j,l)
tw(i,j,l)= - dwdt - (uwx+vwy+wwz) - pz/rhof + (nuf+nuf_rw(i,j,l))*(wxx+wyy+wzz) + Fsz(i,j,l)
! coefficients of u(i,j,l) of momentum eq
Cdwdt= (1._pr+tau01)/(2._pr*tau01*dt)
Cviscz= -2._pr*nuf * ( 1._pr/((xc(i+1)-xc(i))*(xc(i)-xc(i-1))) &
+ 1._pr/((yc(j+1)-yc(j))*(yc(j)-yc(j-1))) &
+ 1._pr/((zf(l+1)-zf(l))*(zf(l)-zf(l-1))) )
Cviscz= -2._pr*(nuf+nuf_rw(i,j,l)) * ( 1._pr/((xc(i+1)-xc(i))*(xc(i)-xc(i-1))) &
+ 1._pr/((yc(j+1)-yc(j))*(yc(j)-yc(j-1))) &
+ 1._pr/((zf(l+1)-zf(l))*(zf(l)-zf(l-1))) )
! coefficients lhs-rhs
qw(i,j,l)= Cdwdt - Cviscz
......@@ -303,3 +301,87 @@
enddo; enddo; enddo
end
!####################################################################
!> @author Holger Grosshans
!> @brief calculate turbulent viscosity with static Smagorinsky model
subroutine turbSmagorinsky
use var
real(kind=pr) :: ux,uy,uz,vx,vy,vz,wx,wy,wz,l_s
integer :: i,j,l
do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
if (celltype(i+1,j,l).ne.wall) then
ux= (u(i+1,j,l)-u(i-1,j,l))/(xf(i+1)-xf(i-1)) ! (2nd order central derivatives)
uy= (u(i,j+1,l)-u(i,j-1,l))/(yc(j+1)-yc(j-1))
uz= (u(i,j,l+1)-u(i,j,l-1))/(zc(l+1)-zc(l-1))
vx= (v(i+1,j,l)+v(i+1,j-1,l)-v(i,j,l)-v(i,j-1,l))/(2._pr*(xc(i+1)-xc(i)))
vy= (v(i+1,j,l)+v(i,j,l)-v(i+1,j-1,l)-v(i,j-1,l))/(2._pr*(yf(j)-yf(j-1)))
vz= (v(i+1,j,l+1)+v(i,j,l+1)+v(i+1,j-1,l+1)+v(i,j-1,l+1)-v(i+1,j,l-1)-v(i,j,l-1)-v(i+1,j-1,l-1)-v(i,j-1,l)-1)/ &
(4._pr*(zc(l+1)-zc(l-1)))
wx= (w(i+1,j,l)+w(i+1,j,l-1)-w(i,j,l)-w(i,j,l-1))/(2._pr*(xc(i+1)-xc(i)))
wy= (w(i+1,j+1,l)+w(i,j+1,l)+w(i+1,j+1,l-1)+w(i,j+1,l-1)-w(i+1,j-1,l)-w(i,j-1,l)-w(i+1,j-1,l-1)-w(i,j-1,l-1))/ &
(4._pr*(yc(j+1)-yc(j-1)))
wz= (w(i+1,j,l)+w(i,j,l)-w(i+1,j,l-1)-w(i,j,l-1))/(2._pr*(zf(l)-zf(l-1)))
l_s= C_s*filteru(i,j,l)
nuf_ru(i,j,l)= l_s**2*S(ux,uy,uz,vx,vy,vz,wx,wy,wz) !< Pope Eq. (13.128)
endif
if (celltype(i,j+1,l).ne.wall) then
ux= (u(i,j+1,l)+u(i,j,l)-u(i-1,j+1,l)-u(i-1,j,l))/(2._pr*(xf(i)-xf(i-1)))
uy= (u(i,j+1,l)+u(i-1,j+1,l)-u(i,j,l)-u(i-1,j,l))/(2._pr*(yc(j+1)-yc(j)))
uz= (u(i-1,j+1,l+1)+u(i,j+1,l+1)+u(i-1,j,l+1)+u(i,j,l+1)-u(i-1,j+1,l-1)-u(i,j+1,l-1)-u(i-1,j,l-1)-u(i,j,l-1))/ &
(4._pr*(zc(l+1)-zc(l-1)))
vx= (v(i+1,j,l)-v(i-1,j,l))/(xc(i+1)-xc(i-1))
vy= (v(i,j+1,l)-v(i,j-1,l))/(yf(j+1)-yf(j-1))
vz= (v(i,j,l+1)-v(i,j,l-1))/(zc(l+1)-zc(l-1))
wx= (w(i+1,j,l)+w(i+1,j+1,l)+w(i+1,j,l-1)+w(i+1,j+1,l-1)-w(i-1,j,l)-w(i-1,j+1,l)-w(i-1,j,l-1)-w(i-1,j+1,l-1))/ &
(4._pr*(xc(i+1)-xc(i-1)))
wy= (w(i,j+1,l)+w(i,j+1,l-1)-w(i,j,l)-w(i,j,l-1))/(2._pr*(yc(j+1)-yc(j)))
wz= (w(i,j,l)+w(i,j+1,l)-w(i,j,l-1)-w(i,j+1,l-1))/(2._pr*(zf(l)-zf(l-1)))
l_s= C_s*filterv(i,j,l)
nuf_rv(i,j,l)= l_s**2*S(ux,uy,uz,vx,vy,vz,wx,wy,wz) !< Pope Eq. (13.128)
endif
if (celltype(i,j,l+1).ne.wall) then
ux= (u(i,j,l)+u(i,j,l+1)-u(i-1,j,l)-u(i-1,j,l+1))/(2._pr*(xf(i)-xf(i-1)))
uy= (u(i,j+1,l)+u(i,j+1,l+1)+u(i-1,j+1,l)+u(i-1,j+1,l+1)-u(i,j-1,l)-u(i,j-1,l+1)-u(i-1,j-1,l)-u(i-1,j-1,l+1))/ &
(4._pr*(yc(j+1)-yc(j-1)))
uz= (u(i,j,l+1)+u(i-1,j,l+1)-u(i,j,l)-u(i-1,j,l))/(2._pr*(zc(l+1)-zc(l)))
vx= (v(i+1,j,l)+v(i+1,j,l+1)+v(i+1,j-1,l)+v(i+1,j-1,l+1)-v(i-1,j,l)-v(i-1,j,l+1)-v(i-1,j-1,l)-v(i-1,j-1,l+1))/ &
(4._pr*(xc(i+1)-xc(i-1)))
vy= (v(i,j,l)+v(i,j,l+1)-v(i,j-1,l)-v(i,j-1,l+1))/(2._pr*(yf(j)-yf(j-1)))
vz= (v(i,j,l+1)+v(i,j-1,l+1)-v(i,j,l)-v(i,j-1,l))/(2._pr*(yc(i+1)-yc(i)))
wx= (w(i+1,j,l)-w(i-1,j,l))/(xc(i+1)-xc(i-1))
wy= (w(i,j+1,l)-w(i,j-1,l))/(yc(j+1)-yc(j-1))
wz= (w(i,j,l+1)-w(i,j,l-1))/(zf(l+1)-zf(l-1))
l_s= C_s*filterw(i,j,l)
nuf_rw(i,j,l)= l_s**2*S(ux,uy,uz,vx,vy,vz,wx,wy,wz) !< Pope Eq. (13.128)
endif
enddo; enddo; enddo
contains
real(kind=pr) function S(ux,uy,uz,vx,vy,vz,wx,wy,wz)
!> @brief calculate characteristic filtered rate of strain, Pope Eq. (13.74)
real(kind=pr) :: ux,uy,uz,vx,vy,vz,wx,wy,wz,Sij(3,3)
Sij(1,1)= ux !< filtered rate of strain tensor, Pope Eq. (2.70)
Sij(2,2)= vy
Sij(3,3)= wz
Sij(1,2)= 0.5_pr*(uy+vx)
Sij(2,1)= Sij(1,2)
Sij(1,3)= 0.5_pr*(uz+wx)
Sij(3,1)= Sij(1,3)
Sij(2,3)= 0.5_pr*(vz+wy)
Sij(3,2)= Sij(2,3)
S= sqrt(2._pr*sum(Sij**2))
end
end
src/main.f90 100644 → 100755
View file @ f40b4e22
File mode changed from 100644 to 100755
src/makefile 100644 → 100755
View file @ f40b4e22
......@@ -9,9 +9,9 @@ INC =
# gcc:
CMP = mpifort
#O3: optimization
#FLAGS = -O3 -mcmodel=medium
FLAGS = -O3 -mcmodel=medium
#debugging:
FLAGS = -g3 -O0 -fbounds-check -mcmodel=medium -fimplicit-none -fcheck=all -fbacktrace -floop-nest-optimize -ffpe-trap=invalid,zero,overflow -Wconversion -fno-tree-vectorize
#FLAGS = -g3 -O0 -fbounds-check -mcmodel=medium -fimplicit-none -fcheck=all -fbacktrace -floop-nest-optimize -ffpe-trap=invalid,zero,overflow -Wconversion -fno-tree-vectorize
#-Wall
FFLAGS = -c $(FLAGS)
......
src/parallel.f90 100644 → 100755
View file @ f40b4e22
......@@ -212,7 +212,7 @@ contains
timenow=real(mpi_wtime(),kind=pr)
! mpi_allreduce does not work with intel compiler
! mpi_allreduce may not work with intel compiler
call mpi_allreduce(myvar,syncMax,1,mpi_pr,mpi_max,mpi_comm_world,mpierr)
timeend=real(mpi_wtime(),kind=pr)
......@@ -220,6 +220,24 @@ contains
end
!####################################################################
!> @author Holger Grosshans
!> @brief compute the min of a scalar over all processors
real(kind=pr) function syncMin(myvar)
use var
use mpi
real(kind=pr) :: myvar
timenow=real(mpi_wtime(),kind=pr)
! mpi_allreduce may not work with intel compiler
call mpi_allreduce(myvar,syncMin,1,mpi_pr,mpi_min,mpi_comm_world,mpierr)
timeend=real(mpi_wtime(),kind=pr)
timecom(1)=timecom(1)+timeend-timenow
end
!####################################################################
!> @author Holger Grosshans
!> @brief compute the sum of a real scalar over all processors
......
src/particles.f90 100644 → 100755
View file @ f40b4e22
!####################################################################
!> @author Holger Grosshans
!> @brief solve particulate phase
subroutine solveParticles
use var
if (pnd.ne.0.or.npTot.ne.0) then
if (((bcx.eq.'i').and.(nt.ge.ntseed)).or.(nt.eq.ntseed)) call particlesSeed
call particlesDrag
call particlesLift
if ((elForceScheme.eq.1).or.(elForceScheme.eq.3)) call forcesGauss
if ((elForceScheme.eq.2).or.(elForceScheme.eq.3)) call forcesCoulomb
call particlesVelocityNext
call particlesCollideNext
call particlesTransportNext
if (myid.eq.0) write(*,'(a,i8,2(a,es8.2e2))') &
'particles |transp. = ',npTot,' |dt/t_el = ',dtNext*rtau_el_max,' |dt/t_p = ',dtNext*rtau_p_max
endif
end
!####################################################################
!> @author Holger Grosshans
!> @brief compute weigths for interpolation from Lagrangian to
......@@ -60,13 +39,13 @@
! volume on the Eulerian grid
if (direction.eq.0) then
volE=(xf(iend)-xf(ibeg-1))*(yf(jend)-yf(jbeg-1))*(zf(lend)-zf(lbeg-1))
volE= sum(vol(ibeg:iend,jbeg:jend,lbeg:lend))
else if (direction.eq.1) then
volE=(xc(iend+1)-xc(ibeg))*(yf(jend)-yf(jbeg-1))*(zf(lend)-zf(lbeg-1))
volE= sum(volu(ibeg:iend,jbeg:jend,lbeg:lend))
else if (direction.eq.2) then
volE=(xf(iend)-xf(ibeg-1))*(yc(jend+1)-yc(jbeg))*(zf(lend)-zf(lbeg-1))
volE= sum(volv(ibeg:iend,jbeg:jend,lbeg:lend))
else if (direction.eq.3) then
volE=(xf(iend)-xf(ibeg-1))*(yf(jend)-yf(jbeg-1))*(zc(lend+1)-zc(lbeg))
volE= sum(volw(ibeg:iend,jbeg:jend,lbeg:lend))
endif
deltaf=xc(iend)-xc(ibeg)
......@@ -100,7 +79,7 @@
npTot=syncSumI(np)
nppTot=syncSumI(npp)
maxnp=npTot*(nrprocs)*3 !if nrprocs=1, one particle can be in domain, gc and sync storage
maxnp=npTot*3 ! one particle can be in domain, gc and sync storage
call allocateParticleArray(nppTot,up)
call allocateParticleArray(nppTot,vp)
......
src/particlesTransport.f90 100644 → 100755
View file @ f40b4e22
!####################################################################
!> @author Holger Grosshans
!> @brief solve particulate phase
subroutine solveParticles
use var
if (pnd.ne.0.or.npTot.ne.0) then
if (((bcx.eq.'i').and.(nt.ge.ntseed)).or.(nt.eq.ntseed)) call particlesSeed
call particlesDrag
call particlesLift
if ((elForceScheme.eq.1).or.(elForceScheme.eq.3)) call forcesGauss
if ((elForceScheme.eq.2).or.(elForceScheme.eq.3)) call forcesCoulomb
call particlesVelocityNext
call particlesCollideNext
call particlesTransportNext
if (myid.eq.0) write(*,'(a,i8,2(a,es8.2e2))') &
'particles |transp. = ',npTot,' |dt/t_el = ',dtNext*rtau_el_max,' |dt/t_p = ',dtNext/tau_p_min
endif
end
!####################################################################
!> @author Holger Grosshans
!> @brief explicit first-order Euler forward integration of particle
......@@ -25,31 +46,31 @@
!####################################################################
!> @author Holger Grosshans
!> @brief particle drag force
!> @brief specific particle drag force
subroutine particlesDrag
use var
use parallel
real(kind=pr) :: rtau_p
real(kind=pr) :: tau_p
integer :: n
call fluidVelocity
rtau_p_max=0._pr
tau_p_min=1.e10_pr
do 1 n=1,np
call calcRtau_p(rtau_p,n)
fx_d(n)= (uf(n)-up(n))*rtau_p
fy_d(n)= (vf(n)-vp(n))*rtau_p
fz_d(n)= (wf(n)-wp(n))*rtau_p
call calcTau_p(tau_p,n)
fx_d(n)= (uf(n)-up(n))/tau_p
fy_d(n)= (vf(n)-vp(n))/tau_p
fz_d(n)= (wf(n)-wp(n))/tau_p
1 enddo
rtau_p_max=syncMax(rtau_p_max)
tau_p_min=syncMin(tau_p_min)
end
!####################################################################
!> @author Holger Grosshans
!> @brief particle lift force (Saffman, 1965, 1968, correction by Mei, 1992)
!> @brief specific particle lift force (Saffman, 1965, 1968, correction by Mei, 1992)
subroutine particlesLift
use var
real(kind=pr) :: urel,Reyp,Reyf,beta,Cls
......@@ -78,25 +99,26 @@
!####################################################################
!> @author Holger Grosshans
!> @brief reciprocal (to avoid NaN) particle response time (Putnam, 1961)
subroutine calcRtau_p(rtau_p,n)
!> @brief particle response time
subroutine calcTau_p(tau_p,n)
use var
real(kind=pr) :: urel,Reyp,Cd,rtau_p
real(kind=pr) :: urel,Reyp,Cd,tau_p
integer :: n
urel= sqrt((uf(n)-up(n))**2+(vf(n)-vp(n))**2+(wf(n)-wp(n))**2)
Reyp= max(1.e-10_pr,2._pr*radp(n)*urel/nuf)
urel= max(1.e-10_pr,urel) ! avoid NaN
Reyp= 2._pr*radp(n)*urel/nuf
if (Reyp.lt.1000._pr) then
if (Reyp.lt.1000._pr) then ! Putnam (1961)
Cd= 24._pr/Reyp * (1._pr + 1._pr/6._pr*Reyp**(2._pr/3._pr))
else
Cd= 0.424_pr
endif
rtau_p= 3._pr*rhof*urel*Cd/(8._pr*rhop*radp(n))
tau_p= 8._pr*rhop*radp(n)/(3._pr*rhof*urel*Cd)
rtau_p_max=max(rtau_p_max,rtau_p)
if (rtau_p.gt.(1._pr/dtNext)) rtau_p= 1._pr/dtNext ! stability criteria Euler forward
tau_p_min=min(tau_p_min,tau_p)
if (tau_p.lt.dtNext) tau_p= dtNext ! stability criteria Euler forward
end
......@@ -399,13 +421,14 @@
npp=np
if (bcx.eq.'i') then
if (myid.eq.0) then
npseed=int(dt*pnd)
npseed=int((t-tstart)*pnd)-npstart
npstart=npstart+npseed
else
npseed=0
endif
else
npseed=int(pnd*dimx*dimy*dimz)
! npseed=2
! npseed=2
endif
npseedTot=syncSumI(npseed)
np=npp+npseed
......
src/post.f90 100644 → 100755
View file @ f40b4e22
......@@ -22,10 +22,10 @@
subroutine monitor
use var
use parallel
real(kind=pr) :: ucl,avu,avv,avw,tke,gradu,graduy,graduz, &
real(kind=pr) :: ucl,avu,avv,avw,gradu,graduy,graduz, &
Rec,Reb,avvp,avwp,avyp,avzp,avvf,avwf,avxp,avqp, &
rmsv,rmsw,C05,A05
integer :: i,j,l,m,N05
integer :: m,N05
logical :: file_ex
character(70) :: filename2
......
src/pre.f90 100644 → 100755
View file @ f40b4e22
......@@ -45,6 +45,8 @@
read(10,*)
read(10,*) int_results,int_restart
read(10,*)
read(10,*) turbModel
read(10,*)
read(10,*) cfl,itmax,tol
read(10,*)
read(10,*) ubulk0,dpdx
......@@ -117,6 +119,8 @@
'Grid in x,y,z-direction [(u)niform/(s)tretch] = ',gridx,gridy,gridz
write(20,'(a,4(i8))') &
'nt: start, compute, results, restart (-,-,-,-) = ',ntstart,ntime,int_results,int_restart
write(20,'(a,3x,a)') &
'turbulence model: = ',turbModel
write(20,'(a,es11.2e2,i11,es11.2e2)') &
'CFL number (-), max number iterations (-), rel. tol. = ',cfl,itmax,tol
write(20,'(a,2(es11.2e2))') &
......@@ -163,11 +167,12 @@
Ex_el(ii,jj,ll),Ey_el(ii,jj,ll),Ez_el(ii,jj,ll), &
u01(ii,jj,ll),v01(ii,jj,ll),w01(ii,jj,ll), &
u02(ii,jj,ll),v02(ii,jj,ll),w02(ii,jj,ll), &
Fsx(ii,jj,ll),Fsy(ii,jj,ll),Fsz(ii,jj,ll))
Fsx(ii,jj,ll),Fsy(ii,jj,ll),Fsz(ii,jj,ll), &
nuf_ru(ii,jj,ll),nuf_rv(ii,jj,ll),nuf_rw(ii,jj,ll))
! particles
allocate(npic(ii,jj,ll))
np=0; npp=0; nseedLoc=0;
np=0; npp=0; nseedLoc=0; npstart=0
call allocateParticleArrays
! init
......@@ -175,9 +180,10 @@
u01=u; v01=v; w01= w
u02=u; v02=v; w02= w
Fsx=0._pr; Fsy=0._pr; Fsz=0._pr
nuf_ru=0._pr; nuf_rv=0._pr; nuf_rw=0._pr
celltype=active
t=0._pr; timecom=0._pr
rtau_p_max=0._pr; Dup_max=0._pr
tau_p_min=0._pr; dup_max=0._pr
vol=0._pr; volu=0._pr; volv=0._pr; volw=0._pr;
Ca=0._pr; Cb1=0._pr; Cb2=0._pr; Cc1=0._pr; Cc2=0._pr; Cd1=0._pr; Cd2=0._pr;
......@@ -217,7 +223,7 @@
subroutine gridGeneration
use var
real(kind=pr) :: hy_temp,hz_temp
integer :: i,j,l,m
integer :: i,j,l
imin=1+gc; imax=ii-gc
jmin=1+gc; jmax=jj-gc
......@@ -335,7 +341,7 @@
use var
integer :: i,j,l
do i=imin,imax; do j=jmin,jmax; do l=lmin,lmax
do i=imin-1,imax+1; do j=jmin-1,jmax+1; do l=lmin-1,lmax+1
vol(i,j,l)= (xf(i)-xf(i-1))*(yf(j)-yf(j-1))*(zf(l)-zf(l-1))
volu(i,j,l)=(xc(i+1)-xc(i))*(yf(j)-yf(j-1))*(zf(l)-zf(l-1))
volv(i,j,l)=(xf(i)-xf(i-1))*(yc(j+1)-yc(j))*(zf(l)-zf(l-1))
......@@ -343,10 +349,15 @@
enddo; enddo; enddo
volTot= sum(vol(imin:imax,jmin:jmax,lmin:lmax))
! LES filter width
filteru= volu**(1._pr/3._pr)
filterv= volv**(1._pr/3._pr)
filterw= volw**(1._pr/3._pr)
! coefficients for PISO algorithm, not sure which version is better
do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
Cb1(i,j,l)=0._pr; Cb2(i,j,l)=0._pr; Cc1(i,j,l)=0._pr
Cc2(i,j,l)=0._pr; Cd1(i,j,l)=0._pr; Cd2(i,j,l)=0._pr
Cb1(i,j,l)=0._pr; Cc1(i,j,l)=0._pr; Cd1(i,j,l)=0._pr
Cb2(i,j,l)=0._pr; Cc2(i,j,l)=0._pr; Cd2(i,j,l)=0._pr
if (celltype(i+1,j,l).ne.wall) Cb1(i,j,l)= 1._pr/((xc(i+1)-xc(i))*(xf(i)-xf(i-1)))
if (celltype(i-1,j,l).ne.wall) Cb2(i,j,l)= 1._pr/((xc(i)-xc(i-1))*(xf(i)-xf(i-1)))
if (celltype(i,j+1,l).ne.wall) Cc1(i,j,l)= 1._pr/((yc(j+1)-yc(j))*(yf(j)-yf(j-1)))
......@@ -355,15 +366,6 @@
if (celltype(i,j,l-1).ne.wall) Cd2(i,j,l)= 1._pr/((zc(l)-zc(l-1))*(zf(l)-zf(l-1)))
Ca(i,j,l)= Cb1(i,j,l)+Cb2(i,j,l)+Cc1(i,j,l)+Cc2(i,j,l)+Cd1(i,j,l)+Cd2(i,j,l)
enddo; enddo; enddo
!do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
! Cb1(i,j,l)= 1._pr/((xc(i+1)-xc(i))*(xf(i)-xf(i-1)))
! Cb2(i,j,l)= 1._pr/((xc(i)-xc(i-1))*(xf(i)-xf(i-1)))
! Cc1(i,j,l)= 1._pr/((yc(j+1)-yc(j))*(yf(j)-yf(j-1)))
! Cc2(i,j,l)= 1._pr/((yc(j)-yc(j-1))*(yf(j)-yf(j-1)))
! Cd1(i,j,l)= 1._pr/((zc(l+1)-zc(l))*(zf(l)-zf(l-1)))
! Cd2(i,j,l)= 1._pr/((zc(l)-zc(l-1))*(zf(l)-zf(l-1)))
! Ca(i,j,l)= Cb1(i,j,l)+Cb2(i,j,l)+Cc1(i,j,l)+Cc2(i,j,l)+Cd1(i,j,l)+Cd2(i,j,l)
!enddo; enddo; enddo
end
......
src/restart.f90 100644 → 100755
View file @ f40b4e22
......@@ -99,5 +99,6 @@
ntstart=ntstart+1
ntend=ntend+1
tstart=t
end
src/timestep.f90 100644 → 100755
View file @ f40b4e22
File mode changed from 100644 to 100755
src/var.f90 100644 → 100755
View file @ f40b4e22
......@@ -20,10 +20,11 @@
real(kind=pr), parameter :: &
eps_el= 8.85e-12_pr, & !< permittivity vacuum/air (F/m)
pi= 4._pr*atan(1._pr), &
urfu= 0.25_pr, & !< under relaxation factor variable u
urfu= 0.25_pr, & !< under-relaxation factor variable u
urfv= urfu, &
urfw= urfv, &
urfp= 4._pr
urfp= 4.0_pr, &
C_s= 0.17_pr !< Smagorinsky constant
real(kind=pr), parameter :: &
Ew=1.e11_pr, & !< duct Young's modulus (kg/s**2/m)
......@@ -66,6 +67,9 @@
dimitot, & !< global number of cells in x-direction
ntseed !< time step particle seeding
character(11) :: &
turbModel !< turbulence model
! run
real(kind=pr) :: &
dpdx, & !< streamwise pressure gradient (N/m**3)
......@@ -74,7 +78,8 @@
dt01, & !< time step size from (nt-2) to (nt-1) (s)
dt02, & !< time step size from (nt-3) to (nt-2) (s)
tau01,tau02, &
t !< physical time (s)
t, & !< physical time (s)
tstart !< start time step
integer :: nt, & !< time step
ntend !< last time step of computation
......@@ -122,11 +127,13 @@
p, & !< fluid pressure (N/m**2)
u01,v01,w01, & !< fluid velocity previous time-step
u02,v02,w02, & !< fluid velocity previous time-step
Fsx,Fsy,Fsz !< source term particles -> fluid (m/s**2)
Fsx,Fsy,Fsz, & !< source term particles -> fluid (m/s**2)
nuf_ru,nuf_rv,nuf_rw, & !< turbulent viscosity (m**2/s)
filteru,filterv,filterw !< LES filter width (m)
! particles
real(kind=pr) :: &
rtau_p_max, & !< max reciprocal particle time-scale in a time-step
tau_p_min, & !< min particle response time in a time-step
rtau_el_max, & !< max reciprocal Coulomb force time-scale in a time-step
dup_max !< max particle velocity change in a time-step
......@@ -135,6 +142,7 @@
npTot, & !< total number of particles
npp, & !< local number of particles before some operation
maxnp, & !< maximum possible local number of particles
npstart, & !< number of particles seeded since tstart
nseedLoc !< local counter for seeded particles
real(kind=pr), allocatable, dimension(:) :: &
......