commit 8c9115e866e9708c1a35e48bac072e9136c861b2 Author: Yeongdo Park Date: Fri Mar 15 10:51:59 2019 -0400 Initial commit diff --git a/code/Compact.f90 b/code/Compact.f90 new file mode 100644 index 0000000..e75fcf1 --- /dev/null +++ b/code/Compact.f90 @@ -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 diff --git a/code/makefile b/code/makefile new file mode 100644 index 0000000..6d1afd7 --- /dev/null +++ b/code/makefile @@ -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 diff --git a/code/test.f90 b/code/test.f90 new file mode 100644 index 0000000..c5b4237 --- /dev/null +++ b/code/test.f90 @@ -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 diff --git a/code/ysolve.f90 b/code/ysolve.f90 new file mode 100644 index 0000000..a96504e --- /dev/null +++ b/code/ysolve.f90 @@ -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