research/dns-hit3d-fdm
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

add scalar equation, disable reaction

Browse source
This commit is contained in:
park0d 2017-03-21 01:03:39 +09:00
parent dac299c060
commit fe6444366c
2 changed files with 95 additions and 61 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
Makefile
View file

@ -4,7 +4,7 @@
MPIF90 = blah
FFTW_HOME = ../libs/fftw-3.2.2
FFTW_HOME = ../libs/fftw
MPIF90 = mpif90
FCFLAGS = -fdefault-double-8 -fdefault-integer-8 -fdefault-real-8 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_HOME)/include
LDFLAGS = -fdefault-double-8 -fdefault-integer-8 -fdefault-real-8 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib -lfftw3 -lm

154
m_fdm_calc.f90
View file

@ -8,6 +8,11 @@ module m_fdm_calc
!variables
integer, parameter :: nscalar = 2
integer, parameter :: neq = 1 + nscalar
integer, parameter :: nd1 = 2 + 2 * nscalar
integer, parameter :: nd2 = nscalar
real*8, dimension(:,:,:), allocatable :: u_,v_,w_
real*8, dimension(:,:,:,:), allocatable :: y1,y2,yf
real*8, dimension(:,:), allocatable :: fz, dfz, fzz, dfzz
@ -234,20 +239,20 @@ module m_fdm_calc
fdmtime=0.
fullsavenum=1000 !full save file
allocate(y1(2,nx,ny,nz))
allocate(y2(2,nx,ny,nz))
allocate(yf(2,nx,ny,nz))
allocate(y1(neq,nx,ny,nz))
allocate(y2(neq,nx,ny,nz))
allocate(yf(neq,nx,ny,nz))
allocate(fz(4*nx*ny,nz))
allocate(dfz(4*nx*ny,nz))
allocate(fzz(nx*ny,nz))
allocate(dfzz(nx*ny,nz))
allocate( fz(nd1*nx*ny,nz))
allocate(dfz(nd1*nx*ny,nz))
allocate( fzz(nd2*nx*ny,nz))
allocate(dfzz(nd2*nx*ny,nz))
allocate(fzu(4*nx*ny,2))
allocate(fzl(4*nx*ny,2))
allocate(fzu(nd1*nx*ny,2))
allocate(fzl(nd1*nx*ny,2))
allocate(fzzu(nx*ny,2))
allocate(fzzl(nx*ny,2))
allocate(fzzu(nd2*nx*ny,2))
allocate(fzzl(nd2*nx*ny,2))
y1=0.0
y2=0.0
@ -280,6 +285,7 @@ module m_fdm_calc
do ii=1,nx
y1(1,ii,i,j)=1. ! rho initializing
y1(2,ii,i,j)=yy(ii) ! Yr initializing
y1(3,ii,i,j)=0.0 ! Yr initializing
enddo
enddo
enddo
@ -467,16 +473,16 @@ module m_fdm_calc
integer :: i,j,k,xx,yy,zz,ii
integer :: n
integer :: idx1, idx2
integer :: idx1, idx2, idx3
real*8 :: wrate,yr,yp
real*8 :: r1_(2,xx,yy,zz),f_(2,xx,yy,zz)
real*8 :: r1_(neq,xx,yy,zz),f_(neq,xx,yy,zz)
real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
real*8 :: ux(4,xx),dux(4,xx),d2ux(xx)
real*8 :: uy(4,yy),duy(4,yy),d2uy(yy)
real*8 :: uz(4,zz),duz(4,zz),d2uz(zz)
real*8 :: uux(xx)
real*8 :: uuy(yy)
real*8 :: uuz(zz)
real*8 :: ux(nd1,xx),dux(nd1,xx),d2ux(nd2,xx)
real*8 :: uy(nd1,yy),duy(nd1,yy),d2uy(nd2,yy)
real*8 :: uz(nd1,zz),duz(nd1,zz),d2uz(nd2,zz)
real*8 :: uux(nd2,xx)
real*8 :: uuy(nd2,yy)
real*8 :: uuz(nd2,zz)
real*8 :: y
@ -485,17 +491,19 @@ module m_fdm_calc
DO j=1,yy
DO i=1,xx
y=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y
! y=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y
wrate=pre*y*exp(-ac/(1.+bc*(1.-y))) !wrate
IF ((1.-y).le.c_ref) THEN
wrate=min_wr
IF ((1.-y).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(1.-y-c_ref))+ &
min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
ENDIF
wrate=0.0
! wrate=pre*y*exp(-ac/(1.+bc*(1.-y))) !wrate
! IF ((1.-y).le.c_ref) THEN
! wrate=min_wr
! IF ((1.-y).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(1.-y-c_ref))+ &
! min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
! ENDIF
f_(1,i,j,k) = 0.0 ! continuity
f_(2,i,j,k) = - wrate ! species conservation
f_(3,i,j,k) = - wrate ! species conservation
ENDDO
ENDDO
@ -509,14 +517,18 @@ module m_fdm_calc
DO i=1,xx
DO k=1,zz
uz(1,k)=r1_(1,i,j,k) ! 1:rho
uz(2,k)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y
uz(2,k)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:YF
uz(3,k)=uz(1,k)*ww_(i,j,k) ! 3:rho*w
uz(4,k)=uz(3,k)*uz(2,k) ! 4:rho*w*Y
uuz (k)=uz(2,k)
uz(4,k)=uz(3,k)*uz(2,k) ! 4:rho*w*YF
uz(5,k)=r1_(3,i,j,k)/r1_(1,i,j,k) ! 5:YO
uz(6,k)=uz(3,k)*uz(5,k) ! 6:rho*w*YO
uuz (1,k)=uz(2,k)
uuz (2,k)=uz(5,k)
ENDDO
CALL dfp(yy,hy,uz(1:4,:),duz(1:4,:),4,3)
CALL d2fp(yy,hy,uuz(:),d2uz(:),1,3)
CALL dfp (yy,hy, uz(1:nd1,:), duz(1:nd1,:),nd1,3)
CALL d2fp(yy,hy,uuz(1:nd2,:),d2uz(1:nd2,:),nd2,3)
DO k=1,zz
! -( d(rho*w)/dz )
@ -525,7 +537,8 @@ module m_fdm_calc
! -( d(rho*w*Yr)/dz ) + d(rho*D* d(Yr)/dz)/dz
! = -( d(rho*w*Yr)/dz )
! + D* (rho* d2(Yr)/dz2 + d(rho)/dz * d(Yr)/dz )
f_(2,i,j,k) = f_(2,i,j,k) - duz(4,k) + diff*(uz(1,k)*d2uz(k)+duz(1,k)*duz(2,k)) ! species conserv.
f_(2,i,j,k) = f_(2,i,j,k) - duz(4,k) + diff*(uz(1,k)*d2uz(1,k)+duz(1,k)*duz(2,k)) ! species conserv.
f_(3,i,j,k) = f_(3,i,j,k) - duz(6,k) + diff*(uz(1,k)*d2uz(2,k)+duz(1,k)*duz(5,k)) ! species conserv.
ENDDO
ENDDO
ENDDO
@ -535,32 +548,39 @@ module m_fdm_calc
DO k=1,zz
DO j=1,yy
DO i=1,xx
idx2 = xx*(j-1)+i
idx1 = (idx2-1)*4
idx3 = xx*(j-1)+i
idx1 = (idx3-1)*nd1
idx2 = (idx3-1)*nd2
fz(idx1+1,k) = r1_(1,i,j,k) ! 1:rho
fz(idx1+2,k) = r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y
fz(idx1+2,k) = r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:YF
fz(idx1+3,k) = r1_(1,i,j,k)*ww_(i,j,k) ! 3:rho*w
fz(idx1+4,k) = r1_(2,i,j,k)*ww_(i,j,k) ! 4:rho*w*Y
fzz(idx2,k) = r1_(2,i,j,k)/r1_(1,i,j,k)
fz(idx1+4,k) = r1_(2,i,j,k)*ww_(i,j,k) ! 4:rho*w*YF
fz(idx1+5,k) = r1_(3,i,j,k)/r1_(1,i,j,k) ! 5:YO
fz(idx1+6,k) = r1_(3,i,j,k)*ww_(i,j,k) ! 6:rho*w*YO
fzz(idx2+1,k) = r1_(2,i,j,k)/r1_(1,i,j,k)
fzz(idx2+2,k) = r1_(3,i,j,k)/r1_(1,i,j,k)
ENDDO
ENDDO
ENDDO
CALL par_dfp (fz, fzl, fzu, dfz, hy, 4*xx*yy, zz, yy, 3)
CALL par_d2fp(fzz, fzzl, fzzu, dfzz, hy, xx*yy, zz, yy, 3)
CALL par_dfp (fz, fzl, fzu, dfz, hy, nd1*xx*yy, zz, yy, 3)
CALL par_d2fp(fzz, fzzl, fzzu, dfzz, hy, nd2*xx*yy, zz, yy, 3)
DO k=1,zz
DO j=1,yy
DO i=1,xx
idx2 = xx*(j-1)+i
idx1 = (idx2-1)*4
idx3 = xx*(j-1)+i
idx1 = (idx3-1)*nd1
idx2 = (idx3-1)*nd2
! -( d(rho*w)/dz )
f_(1,i,j,k) = f_(1,i,j,k) - dfz(idx1+3,k) ! continuity
! -( d(rho*w*Yr)/dz ) + d(rho*D* d(Yr)/dz)/dz
! = -( d(rho*w*Yr)/dz )
! + D* (rho* d2(Yr)/dz2 + d(rho)/dz * d(Yr)/dz )
f_(2,i,j,k) = f_(2,i,j,k) - dfz(idx1+4,k) + diff*(fz(idx1+1,k)*dfzz(idx2,k)+dfz(idx1+1,k)*dfz(idx1+2,k)) ! species conserv.
f_(2,i,j,k) = f_(2,i,j,k) - dfz(idx1+4,k) + diff*(fz(idx1+1,k)*dfzz(idx2+1,k)+dfz(idx1+1,k)*dfz(idx1+2,k)) ! species conserv.
f_(3,i,j,k) = f_(3,i,j,k) - dfz(idx1+6,k) + diff*(fz(idx1+1,k)*dfzz(idx2+2,k)+dfz(idx1+1,k)*dfz(idx1+5,k)) ! species conserv.
ENDDO
ENDDO
@ -573,14 +593,18 @@ module m_fdm_calc
DO j=1,yy
DO i=1,xx
ux(1,i)=r1_(1,i,j,k) ! 1:rho
ux(2,i)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y
ux(2,i)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:YF
ux(3,i)=ux(1,i)*uu_(i,j,k) ! 3:rho*u
ux(4,i)=ux(3,i)*ux(2,i) ! 4:rho*u*Y
uux (i)=ux(2,i)
ux(4,i)=ux(3,i)*ux(2,i) ! 4:rho*u*YF
ux(5,i)=r1_(3,i,j,k)/r1_(1,i,j,k) ! 5:YO
ux(6,i)=ux(3,i)*ux(5,i) ! 6:rho*u*YO
uux (1,i)=ux(2,i)
uux (2,i)=ux(5,i)
ENDDO
CALL dfnonp(xx,hx,ux(1:4,:),dux(1:4,:),4,1)
CALL d2fnonp(xx,hx,uux(:),d2ux(:),1,1)
CALL dfnonp (xx,hx, ux(1:nd1,:), dux(1:nd1,:),nd1,1)
CALL d2fnonp(xx,hx,uux(1:nd2,:),d2ux(1:nd2,:),nd2,1)
DO i=1,xx
@ -590,7 +614,8 @@ module m_fdm_calc
! -( d(rho*u*Yr)/dx ) + d(rho*D* d(Yr)/dx)/dx
! = -( d(rho*u*Yr)/dx )
! + D* (rho* d2(Yr)/dx2 + d(rho)/dx * d(Yr)/dx )
f_(2,i,j,k) = f_(2,i,j,k) - dux(4,i) + diff*(ux(1,i)*d2ux(i)+dux(1,i)*dux(2,i)) ! species conservation
f_(2,i,j,k) = f_(2,i,j,k) - dux(4,i) + diff*(ux(1,i)*d2ux(1,i)+dux(1,i)*dux(2,i)) ! species conservation
f_(3,i,j,k) = f_(3,i,j,k) - dux(6,i) + diff*(ux(1,i)*d2ux(2,i)+dux(1,i)*dux(5,i)) ! species conservation
ENDDO
ENDDO
@ -601,14 +626,18 @@ module m_fdm_calc
DO i=1,xx
DO j=1,yy
uy(1,j)=r1_(1,i,j,k) ! 1:rho
uy(2,j)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:Y
uy(3,j)=uy(1,j)*vv_(i,j,k) ! 3:rho*v
uy(4,j)=uy(3,j)*uy(2,j) ! 4:rho*v*Y
uuy (j)=uy(2,j)
uy(2,j)=r1_(2,i,j,k)/r1_(1,i,j,k) ! 2:YF
uy(3,j)=uy(1,j)*vv_(i,j,k) ! 3:rho*v
uy(4,j)=uy(3,j)*uy(2,j) ! 4:rho*v*YF
uy(5,j)=r1_(3,i,j,k)/r1_(1,i,j,k) ! 5:YO
uy(6,j)=uy(3,j)*uy(5,j) ! 6:rho*v*YO
uuy (1,j)=uy(2,j)
uuy (2,j)=uy(5,j)
ENDDO
CALL dfp(yy,hy,uy(1:4,:),duy(1:4,:),4,2)
CALL d2fp(yy,hy,uuy(:),d2uy(:),1,2)
CALL dfp (yy,hy, uy(1:nd1,:), duy(1:nd1,:),nd1,2)
CALL d2fp(yy,hy,uuy(1:nd2,:),d2uy(1:nd2,:),nd2,2)
DO j=1,yy
! -( d(rho*v)/dy )
@ -617,7 +646,8 @@ module m_fdm_calc
! -( d(rho*v*Yr)/dy ) + d(rho*D* d(Yr)/dy)/dy
! = -( d(rho*v*Yr)/dy )
! + D* (rho* d2(Yr)/dyy2 + d(rho)/dy * d(Yr)/dy )
f_(2,i,j,k)=f_(2,i,j,k)-duy(4,j) + diff*(uy(1,j)*d2uy(j)+duy(1,j)*duy(2,j)) ! species conserv.
f_(2,i,j,k)=f_(2,i,j,k)-duy(4,j) + diff*(uy(1,j)*d2uy(1,j)+duy(1,j)*duy(2,j)) ! species conserv.
f_(3,i,j,k)=f_(3,i,j,k)-duy(6,j) + diff*(uy(1,j)*d2uy(2,j)+duy(1,j)*duy(5,j)) ! species conserv.
ENDDO
ENDDO
ENDDO
@ -628,8 +658,10 @@ module m_fdm_calc
DO j=1,yy
DO i=1,yrsw
f_(2,i,j,k)=r1_(1,i,j,k)*0.+f_(1,i,j,k)*in_yr
f_(3,i,j,k)=r1_(1,i,j,k)*0.+f_(1,i,j,k)*in_yr
ENDDO
IF (uu_(xx,j,k).lt.0.) f_(2,xx,j,k)=f_(1,xx,j,k)*r1_(2,xx,j,k)/r1_(1,xx,j,k)
IF (uu_(xx,j,k).lt.0.) f_(3,xx,j,k)=f_(1,xx,j,k)*r1_(3,xx,j,k)/r1_(1,xx,j,k)
ENDDO
ENDDO
@ -642,7 +674,7 @@ module m_fdm_calc
integer :: istage,xx,yy,zz,i
real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
real*8 :: yy1(2,xx,yy,zz),yy2(2,xx,yy,zz),yyf(2,xx,yy,zz)
real*8 :: yy1(neq,xx,yy,zz),yy2(neq,xx,yy,zz),yyf(neq,xx,yy,zz)
istage=1; CALL substep(yy1,yy1,yy2,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
@ -660,7 +692,7 @@ module m_fdm_calc
integer :: i,j,k,xx,yy,zz
real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
real*8 :: yy1(2,xx,yy,zz),yy2(2,xx,yy,zz),yyf(2,xx,yy,zz)
real*8 :: yy1(neq,xx,yy,zz),yy2(neq,xx,yy,zz),yyf(neq,xx,yy,zz)
CALL RK4(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf)
@ -674,7 +706,7 @@ module m_fdm_calc
integer :: i,j,k,xx,yy,zz,istage
real*8 :: at,bt , wrate , yr
real*8 :: ri(2,xx,yy,zz),r1(2,xx,yy,zz),r2(2,xx,yy,zz),f(2,xx,yy,zz)
real*8 :: ri(neq,xx,yy,zz),r1(neq,xx,yy,zz),r2(neq,xx,yy,zz),f(neq,xx,yy,zz)
real*8 :: a(5),b(5)
real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
integer :: nfinal, iscr, mspec, mpict, msave, nmindt, nv
@ -699,7 +731,7 @@ module m_fdm_calc
DO k=1,zz
DO j=1,yy
DO i=1,xx
DO nv=1,2
DO nv=1,neq
r1(nv,i,j,k)=r1(nv,i,j,k)+at*f(nv,i,j,k)
r2(nv,i,j,k)=r1(nv,i,j,k)+bt*f(nv,i,j,k)
ENDDO
@ -711,15 +743,17 @@ module m_fdm_calc
DO k=1,zz
DO j=1,yy
DO i=1,xx
DO nv=1,2
DO nv=1,neq
r1(nv,i,j,k)=r1(nv,i,j,k)+bt*f(nv,i,j,k)
ENDDO
!==========rho=1 treatment
r1(2,i,j,k)=r1(2,i,j,k)/r1(1,i,j,k)
r1(3,i,j,k)=r1(3,i,j,k)/r1(1,i,j,k)
r1(1,i,j,k)=1.
!==========Max Yr=1 treatment
r1(2,i,j,k)=MIN(in_yr,r1(2,i,j,k))
r1(3,i,j,k)=MIN(in_yr,r1(3,i,j,k))
!==========Min Yr=0 treatment
! r1(2,i,j,k)=MAX(out_yr,r1(2,i,j,k))