Select Git revision
pre.f90 13.51 KiB
!####################################################################
!> @author Holger Grosshans
!> @brief pre processing
subroutine preProcessing
use var
nt=0
call flowConditions
call initVariables
call gridGeneration
call computeCoefficients
call writevtk_grid
if (ntstart.eq.1) then
call initFlow
elseif (ntstart.gt.1) then
call readField
endif
end
!####################################################################
!> @author Holger Grosshans
!> @brief read input data, compute flow conditions and write out
subroutine flowConditions
use var
use mpi
real(kind=pr) :: radpcheck
integer :: dimitotTemp
character :: filename*70,comp_date*8,comp_time*10
! read input file
open(unit=10,file='input/input.dat',status='old')
read(10,*)
read(10,*)
read(10,*)
read(10,*) dimxtot,dimy,dimz
read(10,*)
read(10,*) bcx,bcy,bcz
read(10,*)
read(10,*) dimitotTemp,dimj,diml
read(10,*)
read(10,*) gridx,gridy,gridz
read(10,*)
read(10,*) ntstart,ntime
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
read(10,*)
read(10,*) rhof,nuf
read(10,*)
read(10,*) pnd,ntseed
read(10,*)
read(10,*) radpmin,radpmax,radpstep,rhop,restRatio
read(10,*)
read(10,*) qp0,qpmax
read(10,*)
read(10,*) Qmodel,cSpecies0,Qaccfactor
read(10,*)
read(10,*) g(1),g(2),g(3)
close(10)
radpcheck=mod((radpmax+radpmin)/2._pr,radpstep)
if (myid.eq.0.and.min(abs(radpcheck),abs(radpstep-radpcheck)).lt.1e-10) write(*,'(a,i8)') 'particle radii wrong!'
! connect processors
next=myid+1
prev=myid-1
select case (bcx)
case ('p')
if(prev.eq.-1) prev=nrprocs-1
if(next.eq.nrprocs) next=0
case ('w','i')
if(prev.eq.-1) prev=mpi_proc_null
if(next.eq.nrprocs) next=mpi_proc_null
end select
! compute flow conditions
ntend=ntstart+ntime-1
dimx=dimxtot/nrprocs
dimi=int(dimitotTemp/nrprocs)
dimitot=dimi*nrprocs
if (myid.eq.0.and.dimitot.ne.dimitotTemp) write(*,'(a,i8)') 'dimi set to ',dimitot
npTot=int(pnd*dimxtot*dimy*dimz)
muf=nuf*rhof
ubulk=ubulk0
if (bcy.eq.'w'.and.bcz.eq.'w') delta= min(dimy/2._pr,dimz/2._pr)
if (bcy.eq.'w'.and.bcz.eq.'p') delta= dimy/2._pr
if (bcy.eq.'p'.and.bcz.eq.'w') delta= dimz/2._pr
Re = ubulk*delta/nuf
ftt = dimxtot/ubulk
! write output file
if (myid.ne.0) return
call date_and_time(date=comp_date,time=comp_time)
write(filename,'(3a,a6,a)') 'output/output_',comp_date,'_',comp_time,'.dat'
open(unit=20,file=filename,status='replace')
write(20,'(a,x,a8,a,x,a6)') 'Simulation start date:',comp_date,', time',comp_time
write(20,'(a)') version
write(20,*)
write(20,'(a)') 'Conditions:'
write(20,'(a,i4)') &
'Nr. of processors = ',nrprocs
write(20,'(a,3(es11.2e2))') &
'dimensions in x,y,z-direction (m,m,m) = ',dimxtot,dimy,dimz
write(20,'(a,3x,3(a,x))') &
'BCs in x,y,z-direction [(w)all/(p)eriodic] = ',bcx,bcy,bcz
write(20,'(a,3(i5),a,x,i8)') &
'Nr. of cells in x,y,z-direction (-,-,-) = ',dimitot,dimj,diml,', total =',dimitot*dimj*diml
write(20,'(a,3(i5),a,x,i8)') &
'Nr. of cells in x,y,z-direction per proc. (-,-,-) = ',dimi,dimj,diml,', total =',dimi*dimj*diml
write(20,'(a,3x,3(a,x))') &
'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))') &
'in/initial flow (m/s), pressure gradient (N/m**3) = ',ubulk,dpdx
write(20,'(a,2(es11.2e2))') &
'fluid density (kg/m**3), kinematic visc. (m**2/s) = ',rhof,nuf
write(20,'(a,i11,es11.2,i11)') &
'particle number (-), (-/m**3) or (-/s), seed nt (-) = ',npTot,pnd,ntseed
write(20,'(a,2(es11.2e2))') &
'mat density (kg/m**3), rest coeff (-) = ',rhop,restRatio write(20,'(a,3(es11.2e2))') &
'particle rad min, max, step (m,m,m) = ',radpmin,radpmax,radpstep
write(20,'(a,3(es11.2e2))') &
'particle charge initial, max (C,C), acceleration (-) = ',qp0,qpmax,Qaccfactor
write(20,'(a,3x,a,7x,1(es11.2e2))') &
'charge model, initial species concetration (#/m**2) = ',Qmodel,cSpecies
write(20,'(a,3(es11.2e2))') &
'gravity vector in x,y,z-dir. (m/s**2,m/s**2,m/s**2) = ',g(1),g(2),g(3)
write(20,'(a,3(es11.2e2))') &
'delta (m), Re (-), flow through time (s) = ',delta,Re,ftt
close(20)
end
!####################################################################
!> @brief allocate and initialize all variables
subroutine initVariables
use var
use mpi
ii=dimi+2*gc; jj=dimj+2*gc; ll=diml+2*gc
maxnp=0
! parallel
allocate(mpistatus(mpi_status_size))
! grid
allocate(celltype(ii,jj,ll))
allocate(xc(ii),yc(jj),zc(ll),xf(ii),yf(jj),zf(ll), &
vol(ii,jj,ll),volu(ii,jj,ll),volv(ii,jj,ll),volw(ii,jj,ll), &
Ca(ii,jj,ll),Cb1(ii,jj,ll),Cb2(ii,jj,ll),Cc1(ii,jj,ll),Cc2(ii,jj,ll), &
Cd1(ii,jj,ll),Cd2(ii,jj,ll))
! fluid
allocate(u(ii,jj,ll),v(ii,jj,ll),w(ii,jj,ll),p(ii,jj,ll), &
rho_el(ii,jj,ll),phi_el(ii,jj,ll), &
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), &
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; npstart=0
call allocateParticleArrays
! init
u=0._pr; v=0._pr; w=0._pr; p=0._pr
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
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;
! electrostatics
rho_el=0._pr
phi_el=0._pr
Ex_el= 0._pr; Ey_el=0._pr; Ez_el=0._pr
fx_el= 0._pr; fy_el=0._pr; fz_el=0._pr
fx_d= 0._pr; fy_d= 0._pr; fz_d= 0._pr
fx_l= 0._pr; fy_l= 0._pr; fz_l= 0._pr
q_el= 0._pr; q_elNext=0._pr;
call init_random_seed()
end
!####################################################################
!> @author Holger Grosshans
!> @brief initialize random number generator
subroutine init_random_seed()
use var
integer :: i,n,clock
integer, dimension(:), allocatable :: seed
call random_seed(size=n)
allocate(seed(n))
call system_clock(count=clock)
seed = clock + 37 * (/ (i - 1, i = 1, n) /)
call random_seed(put=seed)
deallocate(seed)
end
!####################################################################
!> @author Holger Grosshans
!> @brief generate grid for gaseous phase
subroutine gridGeneration
use var
real(kind=pr) :: hy_temp,hz_temp
integer :: i,j,l
imin=1+gc; imax=ii-gc
jmin=1+gc; jmax=jj-gc
lmin=1+gc; lmax=ll-gc
gp=ii*jj*ll; dimgp=dimi*dimj*diml; dimgptot=(dimi*nrprocs)*dimj*diml
!> celltypes, default=active
do 1 i=1,ii
select case (bcx)
case ('p')
if ((i.lt.imin).or.(i.gt.imax)) celltype(i,:,:)= passive
case ('w','i')
if ((i.lt.imin).and.(myid.eq.0)) then
celltype(i,:,:)= wall
elseif ((i.lt.imin).and.(myid.ne.0)) then
celltype(i,:,:)= passive
elseif ((i.gt.imax).and.(myid.eq.nrprocs-1)) then
celltype(i,:,:)= wall
elseif ((i.gt.imax).and.(myid.ne.nrprocs-1)) then
celltype(i,:,:)= passive
endif
end select
1 enddo
do 2 j=1,jj
select case (bcy)
case ('p')
if ((j.lt.jmin).or.(j.gt.jmax)) celltype(:,j,:)= passive
case ('w','i')
if ((j.lt.jmin).or.(j.gt.jmax)) celltype(:,j,:)= wall
end select
2 enddo
do 3 l=1,ll
select case (bcz)
case ('p')
if ((l.lt.lmin).or.(l.gt.lmax)) celltype(:,:,l)= passive
case ('w','i')
if ((l.lt.lmin).or.(l.gt.lmax)) celltype(:,:,l)= wall
end select
3 enddo
!> grid
do 4 i=1,ii
if (gridx.eq.'u') then xf(i)=real(i-gc,kind=pr)*dimx/real(dimi,kind=pr)
elseif (gridx.eq.'s') then
!not implemented yet
endif
if (i.eq.1) then
xc(i)= xf(1)-(dimx/real(dimi,kind=pr))/2._pr
else
xc(i)= (xf(i)+xf(i-1))/2._pr
endif
4 enddo
do 5 j=1,jj
if (gridy.eq.'u') then
hy_temp=dimy/real(dimj,kind=pr)
yf(j)=(real(j,kind=pr)-real(jj,kind=pr)/2._pr)*hy_temp
elseif (gridy.eq.'s') then
hy_temp= (1._pr-cos(pi/real(dimj,kind=pr)))*dimy/2._pr !width of first active cell
if (j.lt.jmin) then
yf(j)= real(j+1-jmin,kind=pr)*hy_temp-dimy/2._pr
elseif (j.gt.jmax) then
yf(j)= real(j-jmax,kind=pr)*hy_temp+dimy/2._pr
else
yf(j)=-(cos(real(j-gc,kind=pr)*pi/real(dimj,kind=pr)))*dimy/2._pr
endif
endif
if (j.eq.1) then
yc(j)= yf(1)-hy_temp/2._pr
else
yc(j)= (yf(j)+yf(j-1))/2._pr
endif
5 enddo
do 6 l=1,ll
if (gridz.eq.'u') then
hz_temp=dimz/real(diml,kind=pr)
zf(l)=(real(l,kind=pr)-real(ll,kind=pr)/2._pr)*hz_temp
elseif (gridz.eq.'s') then
hz_temp= (1._pr-cos(pi/real(diml,kind=pr)))*dimz/2._pr !width of first active cell
if (l.lt.lmin) then
zf(l)= real(l+1-lmin,kind=pr)*hz_temp-dimz/2._pr
elseif (l.gt.lmax) then
zf(l)= real(l-lmax,kind=pr)*hz_temp+dimz/2._pr
else
zf(l)=-(cos(real(l-gc,kind=pr)*pi/real(diml,kind=pr)))*dimz/2._pr
endif
endif
if (l.eq.1) then
zc(l)= zf(1)-hz_temp/2._pr
else
zc(l)= (zf(l)+zf(l-1))/2._pr
endif
6 enddo
deltax=xf(ii-gc)*real(myid,kind=pr)
xc=xc+deltax
xf=xf+deltax
xmin=xf(gc); xmax=xf(ii-gc)
ymin=yf(gc); ymax=yf(jj-gc)
zmin=zf(gc); zmax=zf(ll-gc)
end
!####################################################################
!> @author Holger Grosshans
!> @brief pre-compute coefficients to speed of the solution loops subroutine computeCoefficients
use var
integer :: i,j,l
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))
volw(i,j,l)=(xf(i)-xf(i-1))*(yf(j)-yf(j-1))*(zc(l+1)-zc(l))
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; 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)))
if (celltype(i,j-1,l).ne.wall) Cc2(i,j,l)= 1._pr/((yc(j)-yc(j-1))*(yf(j)-yf(j-1)))
if (celltype(i,j,l+1).ne.wall) Cd1(i,j,l)= 1._pr/((zc(l+1)-zc(l))*(zf(l)-zf(l-1)))
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
end
!####################################################################
!> @author Holger Grosshans
!> @brief Initialize the flow to reduce spin-up time (Nelson & Fringer, 2017)
subroutine initFlow
use var
use parallel
real(kind=pr) yplus,random(gp),uplus,vplus,wplus,dist,alpha,utau
integer :: i,j,l,m
call syncRandom(gp,random)
alpha=0.25_pr
m=0
do l=lmin,lmax; do j=jmin,jmax; do i=imin,imax
m=m+1
if (bcy.eq.'w'.and.bcz.eq.'w') dist=delta-max(abs(yc(j)),abs(zc(l)))
if (bcy.eq.'w'.and.bcz.eq.'p') dist=delta-abs(yc(j))
if (bcy.eq.'p'.and.bcz.eq.'w') dist=delta-abs(zc(l))
u(i,j,l)= (1._pr + (2._pr*random(m)-1._pr)*alpha) * ubulk*dist/delta
enddo; enddo; enddo
v=0._pr
w=0._pr
call bcUVW(u,v,w)
u01=u; v01=v; w01=w
u02=u01; v02=v01; w02=w01
end