research/incomp-flame-1d
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Initial commit

Browse source
This commit is contained in:
Yeongdo Park 2019-03-15 10:51:59 -04:00
commit 8c9115e866
4 changed files with 871 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

386
code/Compact.f90 Normal file
View file

@ -0,0 +1,386 @@
MODULE Compact
IMPLICIT NONE
PRIVATE
REAL, DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, &
lyf,lys,wyf,wys, &
lzf,lzs,wzf,wzs
! lyzf,lyzs,wyzf,wyzs
INTEGER :: nxc,nyc,nzc
REAL, PARAMETER :: ezero = 1.0e-14
PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp
CONTAINS
SUBROUTINE ludcmp(nx,ny,nz,xp,yp,zp)
INTEGER, INTENT(IN) :: nx,ny,nz
INTEGER, INTENT(IN) :: xp,yp,zp
INTEGER :: ierr
nxc=nx
nyc=ny
nzc=nz
! IF(nyc /= nzc) PRINT*,'ny should be equal nz'
! xp, yp, zp = 0 : periodic
ALLOCATE(lxf(nxc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
ALLOCATE(lxs(nxc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
IF(xp.eq.0) THEN
ALLOCATE(wxf(nxc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
ALLOCATE(wxs(nxc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
CALL p_lud(1,nxc)
ELSE
CALL nonp_lud(1,nxc)
ENDIF
ALLOCATE(lyf(nyc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
ALLOCATE(lys(nyc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
IF(yp.eq.0) THEN
ALLOCATE(wyf(nyc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
ALLOCATE(wys(nyc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
CALL p_lud(2,nyc)
ELSE
call nonp_lud(2,nyc)
ENDIF
ALLOCATE(lzf(nzc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
ALLOCATE(lzs(nzc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
IF(zp.eq.0) THEN
ALLOCATE(wzf(nzc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
ALLOCATE(wzs(nzc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
CALL p_lud(3,nzc)
ELSE
call nonp_lud(3,nzc)
ENDIF
! CALL x_lud
! CALL yz_lud
END SUBROUTINE ludcmp
SUBROUTINE nonp_lud(xyz,xx)
INTEGER :: i,xyz,xx
REAL, DIMENSION(xx) :: aa
aa=3.
aa(1)=0.5 ; aa(2)=4.
aa(xx-1)=4. ; aa(xx)=0.5
! first derivative
IF (xyz.eq.1) CALL stdlu(aa,xx,lxf) ! x-direction
IF (xyz.eq.2) CALL stdlu(aa,xx,lyf) ! y-direction
IF (xyz.eq.3) CALL stdlu(aa,xx,lzf) ! z-direction
aa=5.5
aa(1)=2./11. ; aa(2)=10.
aa(xx-1)=10. ; aa(xx)=2./11.
! second derivative
IF (xyz.eq.1) CALL stdlu(aa,xx,lxs) ! x-direction
IF (xyz.eq.2) CALL stdlu(aa,xx,lys) ! y-direction
IF (xyz.eq.3) CALL stdlu(aa,xx,lzs) ! z-direction
END SUBROUTINE nonp_lud
SUBROUTINE p_lud(xyz,xx)
INTEGER :: i,xyz,xx
REAL :: a
a=3. ! first derivative
IF (xyz.eq.1) CALL ptdlu(a,xx,lxf,wxf) ! x-direction
IF (xyz.eq.2) CALL ptdlu(a,xx,lyf,wyf) ! y-direction
IF (xyz.eq.3) CALL ptdlu(a,xx,lzf,wzf) ! z-direction
a=11./2. ! second derivative
IF (xyz.eq.1) CALL ptdlu(a,xx,lxs,wxs) ! x-direction
IF (xyz.eq.2) CALL ptdlu(a,xx,lys,wys) ! y-direction
IF (xyz.eq.3) CALL ptdlu(a,xx,lzs,wzs) ! z-direction
END SUBROUTINE p_lud
SUBROUTINE stdlu(a,n,l)
INTEGER :: n
REAL :: a(n),l(n)
REAL :: d
INTEGER :: i
l(1)=1.0/a(1)
DO i=2,n
d=a(i)-l(i-1)
l(i)=1.0/d
ENDDO
END SUBROUTINE stdlu
SUBROUTINE ptdlu(a,n,l,w)
INTEGER :: n
REAL :: a,l(n),w(n)
INTEGER :: i
REAL :: aa(n),d
DO i=1,n-1
aa(i)=a
ENDDO
i=n-1
call stdlu(aa,i,l)
w(1)=1.0
DO i=2,n-2
w(i)=-l(i-1)*w(i-1)
ENDDO
w(n-1)=1.0-l(n-2)*w(n-2)
DO i=1,n-1
w(i)=w(i)*l(i)
ENDDO
d=a
DO i=1,n-1
d=d-w(i)*w(i)/l(i)
ENDDO
l(n)=1./d
END SUBROUTINE ptdlu
SUBROUTINE dfnonp(n,h,x,dx,nd,dir)
INTEGER,INTENT(IN) :: n,nd,dir
REAL,INTENT(IN) :: h
REAL,INTENT(IN),DIMENSION(nd,n) :: x
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j
REAL :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4
h1=1./h
r1=7./3.
r2=1./12.
r3=3.
a=-1.25
b=1.
c=0.25
DO j=1,nd
dx(j,n-1)=x(j,n)-x(j,n-2)
dx(j,n)=-(a*x(j,n)+b*x(j,n-1)+c*x(j,n-2))
dx(j,1)=(a*x(j,1)+b*x(j,2)+c*x(j,3))
dx(j,2)=x(j,3)-x(j,1)
IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) dx(j,n)=0.
IF (x(j,1).eq.x(j,2).and.x(j,2).eq.x(j,3)) dx(j,1)=0.
dx(j,n-1)=dx(j,n-1)*h1*r3
dx(j,n)=dx(j,n)*h1
dx(j,1)=dx(j,1)*h1
dx(j,2)=dx(j,2)*h1*r3
ENDDO
DO i=3,n-2
DO j=1,nd
t1=x(j,i+1)-x(j,i-1)
t2=x(j,i+2)-x(j,i-2)
dx(j,i)=h1*(r1*t1+r2*t2)
ENDDO
ENDDO
IF (dir.eq.1) CALL tdslv(dx,n,lxf,nd) ! x-direction
IF (dir.eq.2) CALL tdslv(dx,n,lyf,nd) ! y-direction
IF (dir.eq.3) CALL tdslv(dx,n,lzf,nd) ! z-direction
END SUBROUTINE dfnonp
SUBROUTINE dfp(n,h,x,dx,nd,dir)
INTEGER,INTENT(IN) :: n,nd,dir
REAL,INTENT(IN) :: h
REAL,INTENT(IN),DIMENSION(nd,n) :: x
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j
REAL :: r1,r2,h1
h1=1./h
r1=7./3.
r2=1./12.
DO j=1,nd
dx(j,n-1)=(r1*(x(j,n)-x(j,n-2))+r2*(x(j,1)-x(j,n-3)))
dx(j,n)=(r1*(x(j,1)-x(j,n-1))+r2*(x(j,2)-x(j,n-2)))
dx(j,1)=(r1*(x(j,2)-x(j,n))+r2*(x(j,3)-x(j,n-1)))
dx(j,2)=(r1*(x(j,3)-x(j,1))+r2*(x(j,4)-x(j,n)))
dx(j,n-1)=dx(j,n-1)*h1
dx(j,n)=dx(j,n)*h1
dx(j,1)=dx(j,1)*h1
dx(j,2)=dx(j,2)*h1
ENDDO
DO i=3,n-2
DO j=1,nd
dx(j,i)=(r1*(x(j,i+1)-x(j,i-1))+r2*(x(j,i+2)-x(j,i-2)))
dx(j,i)=dx(j,i)*h1
ENDDO
ENDDO
IF (dir.eq.1) CALL ptdslv(dx,n,lxf,wxf,nd) ! x-direction
IF (dir.eq.2) CALL ptdslv(dx,n,lyf,wyf,nd) ! y-direction
IF (dir.eq.3) CALL ptdslv(dx,n,lzf,wzf,nd) ! z-direction
END SUBROUTINE dfp
SUBROUTINE ptdslv(r,n,l,w,nd)
INTEGER,INTENT(IN) :: n,nd
REAL,INTENT(INOUT),DIMENSION(nd,n) :: r
REAL,INTENT(IN),DIMENSION(:) :: l,w
INTEGER i,j
REAL, DIMENSION(nd) :: sum
DO j=1,nd
sum(j)=w(1)*r(j,1)
r(j,1)=r(j,1)*l(1)
ENDDO
DO i=2,n-1
DO j=1,nd
r(j,i)=r(j,i)-r(j,i-1)
sum(j)=sum(j)+w(i)*r(j,i)
r(j,i)=r(j,i)*l(i)
ENDDO
ENDDO
DO j=1,nd
r(j,n)=l(n)*(r(j,n)-sum(j))
r(j,n-1)=r(j,n-1)-w(n-1)*r(j,n)
ENDDO
DO i=n-2,1,-1
DO j=1,nd
r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r(j,n)
ENDDO
ENDDO
END SUBROUTINE ptdslv
SUBROUTINE d2fp(n,h,x,dx,nd,dir)
INTEGER,INTENT(IN) :: n,nd,dir
REAL,INTENT(IN) :: h
REAL,INTENT(IN),DIMENSION(nd,n) :: x
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j
REAL :: h2,r1,r2,t1,t2
h2=1./(h*h)
r1=6.
r2=3./8.
DO j=1,nd
t1 = (x(j,n)-2.*x(j,n-1)+x(j,n-2))
t2 = (x(j,1)-2.*x(j,n-1)+x(j,n-3))
IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) t1=0.
IF (x(j,1).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-3)) t2=0.
dx(j,n-1)=(r1*t1+r2*t2)
t1 = (x(j,1)-2.*x(j,n)+x(j,n-1))
t2 = (x(j,2)-2.*x(j,n)+x(j,n-2))
IF (x(j,1).eq.x(j,n).and.x(j,n).eq.x(j,n-1)) t1=0.
IF (x(j,2).eq.x(j,n).and.x(j,n).eq.x(j,n-2)) t2=0.
! dx(j,n)=(r1*(x(j,1)-2.*x(j,n)+x(j,n-1)) &
! +r2*(x(j,2)-2.*x(j,n)+x(j,n-2)))
dx(j,n)=(r1*t1+r2*t2)
t1 = (x(j,2)-2.*x(j,1)+x(j,n))
t2 = (x(j,3)-2.*x(j,1)+x(j,n-1))
IF (x(j,2).eq.x(j,1).and.x(j,1).eq.x(j,n)) t1=0.
IF (x(j,3).eq.x(j,1).and.x(j,1).eq.x(j,n-1)) t2=0.
! dx(j,1)=(r1*(x(j,2)-2.*x(j,1)+x(j,n)) &
! +r2*(x(j,3)-2.*x(j,1)+x(j,n-1)))
dx(j,1)=(r1*t1+r2*t2)
t1 = (x(j,3)-2.*x(j,2)+x(j,1))
t2 = (x(j,4)-2.*x(j,2)+x(j,n))
IF (x(j,3).eq.x(j,2).and.x(j,2).eq.x(j,1)) t1=0.
IF (x(j,4).eq.x(j,2).and.x(j,2).eq.x(j,n)) t2=0.
! dx(j,2)=(r1*(x(j,3)-2.*x(j,2)+x(j,1)) &
! +r2*(x(j,4)-2.*x(j,2)+x(j,n)))
dx(j,2)=(r1*t1+r2*t2)
dx(j,n-1)=dx(j,n-1)*h2
dx(j,n)=dx(j,n)*h2
dx(j,1)=dx(j,1)*h2
dx(j,2)=dx(j,2)*h2
ENDDO
DO i=3,n-2
DO j=1,nd
t1 = (x(j,i+1)-2.*x(j,i)+x(j,i-1))
t2 = (x(j,i+2)-2.*x(j,i)+x(j,i-2))
IF (x(j,i+1).eq.x(j,i).and.x(j,i).eq.x(j,i-1)) t1=0.
IF (x(j,i+2).eq.x(j,i).and.x(j,i).eq.x(j,i-2)) t2=0.
! dx(j,i)=(r1*(x(j,i+1)-2.*x(j,i)+x(j,i-1)) &
! +r2*(x(j,i+2)-2.*x(j,i)+x(j,i-2)))
dx(j,i)=(r1*t1+r2*t2)
dx(j,i)=dx(j,i)*h2
ENDDO
ENDDO
IF (dir.eq.1) CALL ptdslv(dx,n,lxs,wxs,nd) ! x-direction
IF (dir.eq.2) CALL ptdslv(dx,n,lys,wys,nd) ! y-direction
IF (dir.eq.3) CALL ptdslv(dx,n,lzs,wzs,nd) ! z-direction
END SUBROUTINE d2fp
SUBROUTINE tdslv(r,n,l,nd)
INTEGER,INTENT(IN) :: n,nd
REAL,INTENT(INOUT),DIMENSION(nd,n) :: r
REAL,INTENT(IN),DIMENSION(:) :: l
INTEGER i,j
REAL t1
DO j=1,nd
r(j,1)=r(j,1)*l(1)
ENDDO
DO i=2,n
DO j=1,nd
t1=r(j,i)-r(j,i-1)
r(j,i)=l(i)*t1
ENDDO
ENDDO
DO i=n-1,1,-1
DO j=1,nd
r(j,i)=r(j,i)-l(i)*r(j,i+1)
ENDDO
ENDDO
END SUBROUTINE tdslv
SUBROUTINE d2fnonp(n,h,x,dx,nd,dir)
INTEGER,INTENT(IN) :: n,nd,dir
REAL,INTENT(IN) :: h
REAL,INTENT(IN),DIMENSION(nd,n) :: x
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j
REAL :: h2,r1,r2,r3,a,b,c,e,t1,t2
h2=1./(h*h)
r1=6.
r2=3./8.
r3=12.
a=13./11.
b=-27./11.
c=15./11.
e=-1./11.
DO j=1,nd
dx(j,1)=(a*x(j,1)+b*x(j,2)+c*x(j,3)+e*x(j,4))
dx(j,2)=(x(j,3)-2.*x(j,2)+x(j,1))
dx(j,n-1)=(x(j,n)-2.*x(j,n-1)+x(j,n-2))
dx(j,n)=(a*x(j,n)+b*x(j,n-1)+c*x(j,n-2)+e*x(j,n-3))
IF (x(j,1).eq.x(j,2).and.x(j,2).eq.x(j,3).and.x(j,3).eq.x(j,4)) dx(j,1)=0.
IF (x(j,3).eq.x(j,2).and.x(j,2).eq.x(j,1)) dx(j,2)=0.
IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2).and.x(j,n-2).eq.x(j,n-3)) dx(j,n)=0.
IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) dx(j,n-1)=0.
dx(j,1)=dx(j,1)*h2
dx(j,2)=dx(j,2)*h2*r3
dx(j,n-1)=dx(j,n-1)*h2*r3
dx(j,n)=dx(j,n)*h2
ENDDO
DO i=3,n-2
DO j=1,nd
t1 = (x(j,i+1)-2.*x(j,i)+x(j,i-1))
t2 = (x(j,i+2)-2.*x(j,i)+x(j,i-2))
IF (x(j,i+1).eq.x(j,i).and.x(j,i).eq.x(j,i-1)) t1=0.
IF (x(j,i+2).eq.x(j,i).and.x(j,i).eq.x(j,i-2)) t2=0.
! dx(j,i)=(r1*(x(j,i+1)-2.*x(j,i)+x(j,i-1)) &
! +r2*(x(j,i+2)-2.*x(j,i)+x(j,i-2)))
dx(j,i)=(r1*t1+r2*t2)
dx(j,i)=dx(j,i)*h2
ENDDO
ENDDO
IF (dir.eq.1) CALL tdslv(dx,n,lxs,nd) ! x-direction
IF (dir.eq.2) CALL tdslv(dx,n,lys,nd) ! y-direction
IF (dir.eq.3) CALL tdslv(dx,n,lzs,nd) ! z-direction
END SUBROUTINE d2fnonp
END MODULE Compact

17
code/makefile Normal file
View file

@ -0,0 +1,17 @@
flags = -O3 -r8 -i8 -msse3 -axT -m64 -fPIC -i_dynamic
compiler = ifort
ex : test.o ysolve.mod Compact.mod
${compiler} -o ex test.o ysolve.o Compact.o
test.o : test.f90 ysolve.mod
${compiler} -c ${flags} test.f90
ysolve.mod: ysolve.f90 Compact.mod
${compiler} -c ${flags} ysolve.f90
Compact.mod: Compact.f90
${compiler} -c ${flags} Compact.f90
clean:
rm *.o *.mod ex

10
code/test.f90 Normal file
View file

@ -0,0 +1,10 @@
! Premixed Flame Direct Numerical Simulation
! written by Insuk Han
! May/16/2006
! This code is based on the Kim's one
!
PROGRAM test
USE ysolve
IMPLICIT NONE
CALL solve
END PROGRAM test

458
code/ysolve.f90 Normal file
View file

@ -0,0 +1,458 @@
MODULE ysolve
USE Compact
IMPLICIT NONE
PRIVATE
REAL, PARAMETER :: pi=3.14159265358979323846
REAL :: hx,dt,vis,sc,diff,pre,ac,bc,tf,t_now,t_uf,dt_uf
REAL :: c_cut,c_ref,refwr,minf,tar_lo,u0,ctmp,l_0,lo_flm=0.
REAL :: er_lof=0., erdot=0.,min_wr,prof_wr
REAL :: pflame,pflold,oldu
INTEGER :: ncyc=0,int_pr,nx
REAL, DIMENSION(:,:,:), ALLOCATABLE :: u,y1,y2,yf
REAL, DIMENSION(:), ALLOCATABLE :: uxt,duxt
INTEGER :: istage
REAL, DIMENSION(5) :: a=(/ 970286171893./4311952581923., &
6584761158862./12103376702013., &
2251764453980./15575788980749., &
26877169314380./34165994151039., &
0. /), &
b=(/ 1153189308089./22510343858157., &
1772645290293./4653164025191., &
-1672844663538./4480602732383., &
2114624349019./3568978502595., &
5198255086312./14908931495163. /)
PUBLIC :: solve
CONTAINS
!------------------------------------------------------------------------
SUBROUTINE solve
INTEGER :: i,j,k,savenum
REAL :: pflame,pflold,delf=0.
CALL READ_INTRO
CALL ludcmp(nx,1,1,1,0,0)
CALL SET_IC
t_now=0.; t_uf=0.
DO
CALL SET_BC
CALL RK4
IF(t_now.ge.tf) EXIT
IF(t_uf.ge.dt_uf) THEN
pflold=pflame
pflame=0.
DO i=1,nx
pflame=pflame+y1(i,1,1)*hx
uxt(i)=y1(i,1,1)
ENDDO
CALL dfnonp(nx,hx,uxt,duxt,1,1)
delf=1./MAXVAL(ABS(duxt))
oldu=u(1,1,1)
u=u+0.5*(hx*REAL(nx)*tar_lo-pflame)+0.5*(pflold-pflame)
t_uf=0.
WRITE(*,'(a3,f8.3,a10,f6.3,a10,f6.3,a10,f7.4,a10,f7.4,a8,f7.4)') &
' T:',t_now,' // Tar_L:',l_0*tar_lo,' // cur_L:',pflame/hx/REAL(nx)*l_0, &
' // Old_U:',oldu,' // New_U:',u(1,1,1),' // L_f:',delf
! WRITE(*,'(a7,f7.4,a7,f7.4,a10,f7.4)') ' cur_U:',oldu,' // dU:',u(1,1,1)-oldu,' // new_U:',u(1,1,1)
! WRITE(*,*)
ENDIF
ncyc=ncyc+1
t_uf=t_uf+dt
t_now=t_now+dt
IF (MOD(ncyc,int_pr).eq.0) THEN
! WRITE(*,'(a2,f8.3,a9,f10.7,a11,i6,a7,f9.5)') &
! 'T:',t_now,' // dT:',dt,' // NCYC:',ncyc,' // U:',u(1,1,1)
ENDIF
ENDDO
CALL write_sd
CALL write_pre
CALL save_final_field
WRITE(*,*)
WRITE(*,*) 'Fin.'
WRITE(*,*)
END SUBROUTINE solve
SUBROUTINE write_sd
REAL :: c,yr,wrate,dely,sdr,sdd,sd,uu,onelw,dd
INTEGER :: i,j,k,nd
REAL, DIMENSION(1,nx) :: ux,dux,d2ux
REAL, DIMENSION(10,nx) :: sav
sav=0.
! refwr=pre*1.*exp(-ac/(1.+bc*c_cut))
! minf=exp((c_ref-c_cut)*prof_wr)
refwr=pre*1.*exp(-ac/(1.+bc*c_ref))
WRITE(500,*) 'VARIABLES = "X","Yr","C","U","Wrate","|DEL(Y)|","Sdr"'
WRITE(500,*) ' "Sdd","Sd","(1/C)/(dC/dx)","DIV(rho*Dmu*Gra(C))"'
DO i=1,nx
ux(1,i)=y1(i,1,1) ! Yr
IF (ux(1,i).gt.1.) ux(1,i)=1.
ENDDO
nd=1
CALL dfnonp(nx,hx,ux(1,:),dux(1,:),nd,1)
nd=1
CALL d2fnonp(nx,hx,ux(1,:),d2ux(1,:),nd,1)
DO i=1,nx
yr=ux(1,i)
c=1.-yr
IF (c.lt.0.) c=0.
wrate=pre*yr*exp(-ac/(1.+bc*(1.-yr))) !wrate
! IF (c.le.c_cut) THEN
! wrate=min_wr
! IF (c.gt.c_ref) wrate= &
! ((exp((c-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr
! ENDIF
IF (c.le.c_ref) THEN
wrate=min_wr
IF (c.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c-c_ref))+ &
min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
ENDIF
dely=ABS(dux(1,i))
sdr=wrate/dely
sdd=-diff*d2ux(1,i)/dely
IF (dely.eq.0.) THEN
sdr=0.; sdd=0.
ENDIF
sd=sdr+sdd
uu=u(1,1,1)
onelw=(-dux(1,i))/c
dd=-diff*d2ux(1,i)
if (c.eq.0.) onelw=0.
sav(1,i)=sav(1,i)+yr
sav(2,i)=sav(2,i)+c
sav(3,i)=sav(3,i)+uu
sav(4,i)=sav(4,i)+wrate
sav(5,i)=sav(5,i)+dely
sav(6,i)=sav(6,i)+sdr
sav(7,i)=sav(7,i)+sdd
sav(8,i)=sav(8,i)+sd
sav(9,i)=sav(9,i)+onelw
sav(10,i)=sav(10,i)+dd
ENDDO
DO i=1,nx
WRITE(500,'(37e30.20)') (i-1)*hx,sav(1:10,i)
ENDDO
END SUBROUTINE write_sd
SUBROUTINE save_final_field
INTEGER :: i,j,k
OPEN (305,FILE='sfield.bin',form='unformatted',status='unknown')
DO i=1,nx
WRITE (305) y1(i,1,1)
ENDDO
CLOSE (305)
OPEN (305,FILE='sfield.dat')
DO i=1,nx
WRITE (305,'(e30.20)') y1(i,1,1)
ENDDO
CLOSE (305)
END SUBROUTINE save_final_field
SUBROUTINE write_pre
REAL :: yr,c,dy,maxdy=0.,del_f
REAL :: S_L=0.,wrate
INTEGER :: i
REAL, DIMENSION(1,nx) :: ux, dux
DO i=1,nx
yr=y1(i,1,1)
c=1.-yr
ux(1,i)=yr
wrate=pre*yr*exp(-ac/(1.+bc*(1.-yr))) !wrate
! IF (c.le.c_cut) THEN
! wrate=min_wr
! IF (c.gt.c_ref) wrate= &
! ((exp((c-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr
! ENDIF
IF (c.le.c_ref) THEN
wrate=min_wr
IF (c.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c-c_ref))+ &
min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
ENDIF
S_L=S_L+wrate*hx
ENDDO
WRITE(*,'(a31,e14.8)') ' INTEGRAL( Wrate x dx ) => Sc :',S_L
CALL dfnonp(nx,hx,ux(1,:),dux(1,:),1,1)
DO i=(nx/10),nx-(nx/10)
dy=dux(1,i)
maxdy=MAX(ABS(dy),maxdy)
ENDDO
del_f=1./maxdy
WRITE(*,'(a13,e14.8,a25,e14.8)') ' Grid size : ',hx,' / Laminar flame speed : ',u(1,1,1)
WRITE(*,'(a19,e14.8,a3,f9.5,a20)') &
' Flame thickness : ',del_f,' / ',del_f/hx,' grids in the flame.'
WRITE(*,*)
END SUBROUTINE write_pre
SUBROUTINE SET_BC
y1(1,1,1)=1.-ctmp
END SUBROUTINE SET_BC
SUBROUTINE SET_IC
INTEGER :: i, ifl, si
REAL :: xi
ifl=INT(2.0*pi*0.1/hx)
u=u0; si=INT(nx*(1.-tar_lo))
DO i=1,nx
IF(i< nx-(si+ifl/2)) THEN
xi=0.+ctmp
ELSE IF(i> nx-(si-ifl/2)) THEN
xi=1.
ELSE
xi=0.5+REAL(i-nx+si)/REAL(ifl)
ENDIF
! y1(i,1,1)=(1.-xi)*1.
y1(i,1,1)=(1.-xi) ! reactant mass fraction
ENDDO
pflame=0.
DO i=1,nx
pflame=pflame+y1(i,1,1)*hx
ENDDO
END SUBROUTINE SET_IC
SUBROUTINE READ_INTRO
CHARACTER(LEN=8) :: cdum
INTEGER :: itape=300,otape=301,ierr
OPEN(itape,FILE='itape')
OPEN(otape,FILE='otape')
READ(itape,*) cdum,nx
WRITE(otape,*) cdum,nx
READ(itape,*) cdum,l_0
WRITE(otape,*) cdum,l_0
READ(itape,*) cdum,int_pr
WRITE(otape,*) cdum,int_pr
READ(itape,*) cdum,tar_lo
WRITE(otape,*) cdum,tar_lo
READ(itape,*) cdum,dt
WRITE(otape,*) cdum,dt
READ(itape,*) cdum,sc
WRITE(otape,*) cdum,sc
READ(itape,*) cdum,vis
WRITE(otape,*) cdum,vis
READ(itape,*) cdum,pre
WRITE(otape,*) cdum,pre
READ(itape,*) cdum,ac
WRITE(otape,*) cdum,ac
READ(itape,*) cdum,bc
WRITE(otape,*) cdum,bc
READ(itape,*) cdum,u0
WRITE(otape,*) cdum,u0
READ(itape,*) cdum,tf
WRITE(otape,*) cdum,tf
READ(itape,*) cdum,dt_uf
WRITE(otape,*) cdum,dt_uf
READ(itape,*) cdum,ctmp
WRITE(otape,*) cdum,ctmp
READ(itape,*) cdum,c_cut
WRITE(otape,*) cdum,c_cut
! READ(itape,*) cdum,cs
! WRITE(otape,*) cdum,cs
READ(itape,*) cdum,c_ref
WRITE(otape,*) cdum,c_ref
READ(itape,*) cdum,min_wr
WRITE(otape,*) cdum,min_wr
READ(itape,*) cdum,prof_wr
WRITE(otape,*) cdum,prof_wr
CLOSE(itape)
! hx=l_0*pi/REAL(nx)
hx=l_0*pi/REAL(nx-1)
cdum='hx'
WRITE(otape,*) cdum,hx
cdum='Ta/Tu'
WRITE(otape,*) cdum,ac
cdum='Tb/Tu'
WRITE(otape,*) cdum,bc+1
diff=vis/sc
cdum='diff'
WRITE(otape,*) cdum,diff
! refwr=pre*1.*exp(-ac/(1.+bc*c_cut))
! minf=exp((c_ref-c_cut)*prof_wr)
refwr=pre*1.*exp(-ac/(1.+bc*c_ref))
l_0=l_0*pi
ALLOCATE(u(nx,1,1),STAT=ierr) ; u=0.
ALLOCATE(y1(nx,1,1),STAT=ierr) ; y1=0.
ALLOCATE(y2(nx,1,1),STAT=ierr) ; y2=0.
ALLOCATE(yf(nx,1,1),STAT=ierr) ; yf=0.
ALLOCATE(uxt(nx),STAT=ierr) ; uxt=0.
ALLOCATE(duxt(nx),STAT=ierr) ; duxt=0.
END SUBROUTINE READ_INTRO
SUBROUTINE RK4
istage=1; CALL substep(y1,y1,y2,yf)
istage=2; CALL substep(y1,y2,y1,yf)
istage=3; CALL substep(y2,y1,y2,yf)
istage=4; CALL substep(y1,y2,y1,yf)
istage=5; CALL substep(y2,y1,y2,yf)
END SUBROUTINE RK4
!------------------------------------------------------------------------
SUBROUTINE substep(ri,r1,r2,f)
REAL, INTENT(INOUT),DIMENSION(:,:,:) :: ri,r1,r2
REAL, INTENT(OUT),DIMENSION(:,:,:) :: f
INTEGER :: i,j,k
REAL :: at,bt
CALL fns(ri,f)
IF(istage<5) THEN
at=a(istage)*dt
bt=(b(istage)-a(istage))*dt
DO k=1,1 ! nz
DO j=1,1 ! ny
DO i=1,nx
r1(i,j,k)=r1(i,j,k)+at*f(i,j,k)
r2(i,j,k)=r1(i,j,k)+bt*f(i,j,k)
ENDDO
ENDDO
ENDDO
ELSE
bt=b(istage)*dt
DO k=1,1 ! nz
DO j=1,1 ! ny
DO i=1,nx
r1(i,j,k)=r1(i,j,k)+bt*f(i,j,k)
ENDDO
ENDDO
ENDDO
ENDIF
END SUBROUTINE substep
!------------------------------------------------------------------------
SUBROUTINE fns(r1,f)
REAL, INTENT(IN),DIMENSION(:,:,:) :: r1
REAL, INTENT(OUT),DIMENSION(:,:,:) :: f
REAL, DIMENSION(3,nx) :: ux,dux,d2ux
INTEGER :: i,j,k
REAL :: wrate,Ly,Dy
DO k=1,1 !nz
! x-direction
DO j=1,1 !ny
DO i=1,nx
ux(1,i)=u(i,j,k)*r1(i,j,k) ! u*Y
ux(2,i)=u(i,j,k) ! u
ux(3,i)=r1(i,j,k) ! Y
ENDDO
CALL dfnonp(nx,hx,ux(1:3,:),dux(1:3,:),3,1)
CALL d2fnonp(nx,hx,ux(3:3,:),d2ux(1,:),1,1)
DO i=1,nx
wrate=pre*ux(3,i)*exp(-ac/(1.+bc*(1.-ux(3,i)))) !wrate
! IF ((1.-ux(3,i)).le.c_cut) THEN
! wrate=min_wr
! IF ((1.-ux(3,i)).gt.c_ref) wrate= &
! ((exp(((1.-ux(3,i))-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr
! ENDIF
IF ((1.-ux(3,i)).le.c_ref) THEN
wrate=min_wr
IF ((1.-ux(3,i)).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(1.-ux(3,i)-c_ref))+ &
min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
ENDIF
! -0.5*( d(u*Y)/dx + u*dY/dx + Y*du/dx ) + D*d2Y/d2x
f(i,j,k)=-0.5*( dux(1,i) + ux(2,i)*dux(3,i) + &
ux(3,i)*dux(2,i) ) &
+ diff*d2ux(1,i) - wrate
IF (i.eq.nx) THEN
Ly=ux(2,nx)*dux(3,nx) ! Ly = u*dYr/dx
f(nx,1,1)=-wrate
ENDIF
ENDDO
ENDDO
!! y-direction
! DO i=1,nx
! DO j=1,ny
! uy(1,j)=v(i,j,k)*r1(i,j,k) ! v*Y
! uy(2,j)=v(i,j,k) ! v
! uy(3,j)=r1(i,j,k) ! Y
! ENDDO
! CALL dfyz3(ny,hy,uy(1:3,:),duy)
! CALL d2fyz1(ny,hy,uy(3:3,:),d2uy)
! DO j=1,ny
! f(i,j,k)=f(i,j,k) + &
! -0.5*( duy(1,j) + uy(2,j)*duy(3,j) + uy(3,j)*duy(2,j) ) &
! + diff*d2uy(j)
! ENDDO
! ENDDO
ENDDO
! f(1,1,1)=1.0
! f(nx,1,1)=-0.5*( dux(1,nx) + ux(2,nx)*dux(3,nx) + ux(3,nx)*dux(2,nx) )
! f(nx,1,1)=0.0
!! z-direction
! DO j=1,ny
! DO i=1,nx
! DO k=1,nz
! uy(1,k)=w(i,j,k)*r1(i,j,k) ! w*Y
! uy(2,k)=w(i,j,k) ! w
! uy(3,k)=r1(i,j,k) ! Y
! ENDDO
! CALL dfyz3(ny,hy,uy(1:3,:),duy)
! CALL d2fyz1(ny,hy,uy(3:3,:),d2uy)
! DO k=1,nz
! f(i,j,k)=f(i,j,k) + &
! -0.5*( duy(1,k) + uy(2,k)*duy(3,k) + uy(3,k)*duy(2,k) ) &
! + diff*d2uy(k)
! ENDDO
! ENDDO
! ENDDO
! DO j=1,ny
! DO k=1,nz
! f(1,j,k)=0.0
! f(nx,j,k)=-0.5*(dux(1,nx)+ux(2,nx)*dux(3,nx)+ux(3,nx)*&
! dux(2,nx))
! ENDDO
! ENDDO
! Boundary conditionS
f(1,1,1)=0.
Dy=Ly
f(nx,1,1)=f(nx,1,1)-Dy
END SUBROUTINE fns
!------------------------------------------------------------------------
END MODULE ysolve