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

Compare commits

base: research:two-step
Branches Tags
research:master
research:two-step
research:sutherland
research:nonunityle
research:stdin
research:variable_d
..
compare: research:master
Branches Tags
research:two-step
research:sutherland
research:nonunityle
research:stdin
research:variable_d
research:master

No commits in common. "two-step" and "master" have entirely different histories.
two-step ... master

12 changed files with 355 additions and 1202 deletions
Show stats Expand all files Collapse all files

4
.gitlab-ci.yml
View file

@ -16,6 +16,4 @@ test:
script: script:
- cd ${CI_PROJECT_DIR} - cd ${CI_PROJECT_DIR}
- cd code - cd code
- make clean && make - make test
- cd sample/4pi-IC1
- ../../ex

259
code/m_compact.f90 → code/Compact.f90
View file

@ -1,13 +1,14 @@
MODULE m_compact MODULE Compact
IMPLICIT NONE IMPLICIT NONE
PRIVATE
REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, & REAL, DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, &
lyf,lys,wyf,wys, & lyf,lys,wyf,wys, &
lzf,lzs,wzf,wzs lzf,lzs,wzf,wzs
! lyzf,lyzs,wyzf,wyzs
INTEGER :: nxc,nyc,nzc INTEGER :: nxc,nyc,nzc
REAL, PARAMETER :: ezero = 1.0e-14
REAL(KIND=8), PARAMETER :: ezero = 1.0e-14 PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp
CONTAINS CONTAINS
@ -19,22 +20,6 @@
nxc=nx nxc=nx
nyc=ny nyc=ny
nzc=nz nzc=nz
CALL ludcmp_allocate(nx,ny,nz,xp,yp,zp)
CALL ludcmp_calculate(nx,ny,nz,xp,yp,zp)
END SUBROUTINE ludcmp
SUBROUTINE ludcmp_allocate(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' ! IF(nyc /= nzc) PRINT*,'ny should be equal nz'
! xp, yp, zp = 0 : periodic ! xp, yp, zp = 0 : periodic
@ -47,6 +32,9 @@
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
ALLOCATE(wxs(nxc),STAT=ierr) ALLOCATE(wxs(nxc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
CALL p_lud(1,nxc)
ELSE
CALL nonp_lud(1,nxc)
ENDIF ENDIF
ALLOCATE(lyf(nyc),STAT=ierr) ALLOCATE(lyf(nyc),STAT=ierr)
@ -58,6 +46,9 @@
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
ALLOCATE(wys(nyc),STAT=ierr) ALLOCATE(wys(nyc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
CALL p_lud(2,nyc)
ELSE
call nonp_lud(2,nyc)
ENDIF ENDIF
ALLOCATE(lzf(nzc),STAT=ierr) ALLOCATE(lzf(nzc),STAT=ierr)
@ -69,140 +60,19 @@
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
ALLOCATE(wzs(nzc),STAT=ierr) ALLOCATE(wzs(nzc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
ENDIF
END SUBROUTINE ludcmp_allocate
SUBROUTINE ludcmp_deallocate(xp,yp,zp)
INTEGER, INTENT(IN) :: xp,yp,zp
! IF(nyc /= nzc) PRINT*,'ny should be equal nz'
! xp, yp, zp = 0 : periodic
DEALLOCATE(lxf)
DEALLOCATE(lxs)
IF(xp.eq.0) THEN
DEALLOCATE(wxf)
DEALLOCATE(wxs)
ENDIF
DEALLOCATE(lyf)
DEALLOCATE(lys)
IF(yp.eq.0) THEN
DEALLOCATE(wyf)
DEALLOCATE(wys)
ENDIF
DEALLOCATE(lzf)
DEALLOCATE(lzs)
IF(zp.eq.0) THEN
DEALLOCATE(wzf)
DEALLOCATE(wzs)
ENDIF
END SUBROUTINE ludcmp_deallocate
SUBROUTINE ludcmp_testalloc
IF (.not. ALLOCATED(lxf)) print *, "lxf not allocated"
IF (.not. ALLOCATED(lxs)) print *, "lxs not allocated"
IF (.not. ALLOCATED(wxf)) print *, "wxf not allocated"
IF (.not. ALLOCATED(wxs)) print *, "wxs not allocated"
IF (.not. ALLOCATED(lyf)) print *, "lyf not allocated"
IF (.not. ALLOCATED(lys)) print *, "lys not allocated"
IF (.not. ALLOCATED(wyf)) print *, "wyf not allocated"
IF (.not. ALLOCATED(wys)) print *, "wys not allocated"
IF (.not. ALLOCATED(lzf)) print *, "lzf not allocated"
IF (.not. ALLOCATED(lzs)) print *, "lzs not allocated"
IF (.not. ALLOCATED(wzf)) print *, "wzf not allocated"
IF (.not. ALLOCATED(wzs)) print *, "wzs not allocated"
END SUBROUTINE ludcmp_testalloc
SUBROUTINE ludcmp_calculate(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
! CALL ludcmp_testalloc
! IF(nyc /= nzc) PRINT*,'ny should be equal nz'
! xp, yp, zp = 0 : periodic
IF(xp.eq.0) THEN
CALL p_lud(1,nxc)
ELSE
CALL nonp_lud(1,nxc)
ENDIF
IF(yp.eq.0) THEN
CALL p_lud(2,nyc)
ELSE
call nonp_lud(2,nyc)
ENDIF
IF(zp.eq.0) THEN
CALL p_lud(3,nzc) CALL p_lud(3,nzc)
ELSE ELSE
call nonp_lud(3,nzc) call nonp_lud(3,nzc)
ENDIF ENDIF
END SUBROUTINE ludcmp_calculate ! CALL x_lud
! CALL yz_lud
SUBROUTINE test_nonp_lud1(xx, coef)
INTEGER :: xx
REAL(KIND=8), DIMENSION(xx) :: aa
REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef
aa=3.
aa(1)=0.5 ; aa(2)=4.
aa(xx-1)=4. ; aa(xx)=0.5
CALL stdlu(aa,xx,coef)
END SUBROUTINE test_nonp_lud1
SUBROUTINE test_nonp_lud2(xx, coef)
INTEGER :: xx
REAL(KIND=8), DIMENSION(xx) :: aa
REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef
aa=5.5
aa(1)=2./11. ; aa(2)=10.
aa(xx-1)=10. ; aa(xx)=2./11.
CALL stdlu(aa,xx,coef)
END SUBROUTINE test_nonp_lud2
SUBROUTINE test_p_lud1(xx, coef1, coef2)
INTEGER :: xx
REAL(KIND=8) :: a
REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef1, coef2
a=3. ! first derivative
CALL ptdlu(a,xx,coef1,coef2) ! x-direction
END SUBROUTINE test_p_lud1
SUBROUTINE test_p_lud2(xx, coef1, coef2)
INTEGER :: xx
REAL(KIND=8) :: a
REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef1, coef2
a=11./2. ! second derivative
CALL ptdlu(a,xx,coef1,coef2) ! x-direction
END SUBROUTINE test_p_lud2
END SUBROUTINE ludcmp
SUBROUTINE nonp_lud(xyz,xx) SUBROUTINE nonp_lud(xyz,xx)
INTEGER :: i,xyz,xx INTEGER :: i,xyz,xx
REAL(KIND=8), DIMENSION(xx) :: aa REAL, DIMENSION(xx) :: aa
aa=3. aa=3.
aa(1)=0.5 ; aa(2)=4. aa(1)=0.5 ; aa(2)=4.
aa(xx-1)=4. ; aa(xx)=0.5 aa(xx-1)=4. ; aa(xx)=0.5
@ -222,7 +92,7 @@
SUBROUTINE p_lud(xyz,xx) SUBROUTINE p_lud(xyz,xx)
INTEGER :: i,xyz,xx INTEGER :: i,xyz,xx
REAL(KIND=8) :: a REAL :: a
a=3. ! first derivative a=3. ! first derivative
IF (xyz.eq.1) CALL ptdlu(a,xx,lxf,wxf) ! x-direction 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.2) CALL ptdlu(a,xx,lyf,wyf) ! y-direction
@ -235,23 +105,21 @@
SUBROUTINE stdlu(a,n,l) SUBROUTINE stdlu(a,n,l)
INTEGER :: n INTEGER :: n
REAL(KIND=8), INTENT(IN) :: a(n) REAL :: a(n),l(n)
REAL(KIND=8), INTENT(OUT) :: l(n) REAL :: d
REAL(KIND=8) :: d
INTEGER :: i INTEGER :: i
l(1)=1.0d0/a(1) l(1)=1.0/a(1)
DO i=2,n DO i=2,n
d=a(i)-l(i-1) d=a(i)-l(i-1)
l(i)=1.0d0/d l(i)=1.0/d
ENDDO ENDDO
END SUBROUTINE stdlu END SUBROUTINE stdlu
SUBROUTINE ptdlu(a,n,l,w) SUBROUTINE ptdlu(a,n,l,w)
INTEGER :: n INTEGER :: n
REAL(KIND=8), INTENT(OUT) :: a REAL :: a,l(n),w(n)
REAL(KIND=8), INTENT(OUT) :: l(n),w(n)
INTEGER :: i INTEGER :: i
REAL(KIND=8) :: aa(n),d REAL :: aa(n),d
DO i=1,n-1 DO i=1,n-1
aa(i)=a aa(i)=a
@ -273,19 +141,18 @@
l(n)=1./d l(n)=1./d
END SUBROUTINE ptdlu END SUBROUTINE ptdlu
SUBROUTINE dfnonp(n,h,x,dx,nd,dir)
SUBROUTINE rhs1np(n,h,x,dx,nd) INTEGER,INTENT(IN) :: n,nd,dir
INTEGER,INTENT(IN) :: n,nd REAL,INTENT(IN) :: h
REAL(KIND=8),INTENT(IN) :: h REAL,INTENT(IN),DIMENSION(nd,n) :: x
REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j INTEGER :: i,j
REAL(KIND=8) :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4 REAL :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4
h1=1.d0/h h1=1./h
r1=7.d0/3.d0 r1=7./3.
r2=1.d0/12.d0 r2=1./12.
r3=3. r3=3.
a=-1.25 a=-1.25
b=1. b=1.
@ -311,20 +178,6 @@
dx(j,i)=h1*(r1*t1+r2*t2) dx(j,i)=h1*(r1*t1+r2*t2)
ENDDO ENDDO
ENDDO ENDDO
END SUBROUTINE rhs1np
SUBROUTINE dfnonp(n,h,x,dx,nd,dir)
INTEGER,INTENT(IN) :: n,nd,dir
REAL(KIND=8),INTENT(IN) :: h
REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j
REAL(KIND=8) :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4
CALL rhs1np (n,h,x,dx,nd)
IF (dir.eq.1) CALL tdslv(dx,n,lxf,nd) ! x-direction 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.2) CALL tdslv(dx,n,lyf,nd) ! y-direction
IF (dir.eq.3) CALL tdslv(dx,n,lzf,nd) ! z-direction IF (dir.eq.3) CALL tdslv(dx,n,lzf,nd) ! z-direction
@ -332,14 +185,11 @@
SUBROUTINE dfp(n,h,x,dx,nd,dir) SUBROUTINE dfp(n,h,x,dx,nd,dir)
INTEGER,INTENT(IN) :: n,nd,dir INTEGER,INTENT(IN) :: n,nd,dir
REAL(KIND=8),INTENT(IN) :: h REAL,INTENT(IN) :: h
REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x REAL,INTENT(IN),DIMENSION(nd,n) :: x
REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j INTEGER :: i,j
REAL(KIND=8) :: r1,r2,h1 REAL :: r1,r2,h1
! print *, "dfnonp received (nd,n)", nd, n
h1=1./h h1=1./h
r1=7./3. r1=7./3.
@ -371,10 +221,10 @@
SUBROUTINE ptdslv(r,n,l,w,nd) SUBROUTINE ptdslv(r,n,l,w,nd)
INTEGER,INTENT(IN) :: n,nd INTEGER,INTENT(IN) :: n,nd
REAL(KIND=8),INTENT(INOUT),DIMENSION(nd,n) :: r REAL,INTENT(INOUT),DIMENSION(nd,n) :: r
REAL(KIND=8),INTENT(IN),DIMENSION(n) :: l,w REAL,INTENT(IN),DIMENSION(:) :: l,w
INTEGER i,j INTEGER i,j
REAL(KIND=8), DIMENSION(nd) :: sum REAL, DIMENSION(nd) :: sum
DO j=1,nd DO j=1,nd
sum(j)=w(1)*r(j,1) sum(j)=w(1)*r(j,1)
r(j,1)=r(j,1)*l(1) r(j,1)=r(j,1)*l(1)
@ -399,13 +249,11 @@
SUBROUTINE d2fp(n,h,x,dx,nd,dir) SUBROUTINE d2fp(n,h,x,dx,nd,dir)
INTEGER,INTENT(IN) :: n,nd,dir INTEGER,INTENT(IN) :: n,nd,dir
REAL(KIND=8),INTENT(IN) :: h REAL,INTENT(IN) :: h
REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x REAL,INTENT(IN),DIMENSION(nd,n) :: x
REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j INTEGER :: i,j
REAL(KIND=8) :: h2,r1,r2,t1,t2 REAL :: h2,r1,r2,t1,t2
h2=1./(h*h) h2=1./(h*h)
r1=6. r1=6.
r2=3./8. r2=3./8.
@ -464,10 +312,10 @@
SUBROUTINE tdslv(r,n,l,nd) SUBROUTINE tdslv(r,n,l,nd)
INTEGER,INTENT(IN) :: n,nd INTEGER,INTENT(IN) :: n,nd
REAL(KIND=8),INTENT(INOUT),DIMENSION(nd,n) :: r REAL,INTENT(INOUT),DIMENSION(nd,n) :: r
REAL(KIND=8),INTENT(IN),DIMENSION(n) :: l REAL,INTENT(IN),DIMENSION(:) :: l
INTEGER i,j INTEGER i,j
REAL(KIND=8) t1 REAL t1
DO j=1,nd DO j=1,nd
r(j,1)=r(j,1)*l(1) r(j,1)=r(j,1)*l(1)
ENDDO ENDDO
@ -486,12 +334,11 @@
SUBROUTINE d2fnonp(n,h,x,dx,nd,dir) SUBROUTINE d2fnonp(n,h,x,dx,nd,dir)
INTEGER,INTENT(IN) :: n,nd,dir INTEGER,INTENT(IN) :: n,nd,dir
REAL(KIND=8),INTENT(IN) :: h REAL,INTENT(IN) :: h
REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x REAL,INTENT(IN),DIMENSION(nd,n) :: x
REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j INTEGER :: i,j
REAL(KIND=8) :: h2,r1,r2,r3,a,b,c,e,t1,t2 REAL :: h2,r1,r2,r3,a,b,c,e,t1,t2
h2=1./(h*h) h2=1./(h*h)
r1=6. r1=6.
@ -536,4 +383,4 @@
END SUBROUTINE d2fnonp END SUBROUTINE d2fnonp
END MODULE m_compact END MODULE Compact

154
code/m_chemistry.f90
View file

@ -1,154 +0,0 @@
module m_chemistry
use m_parameters
implicit none
real, private :: coef(10)
real, private :: lambda_onestep
real, private :: lambda1_twostep
real, private :: lambda2_twostep
real, private :: beta1_twostep
real, private :: hrp_twostep
contains
subroutine init_chemistry
character(len=40) :: nrxn_string
if (nrxn == 1) then
reaction_type = "onestep"
else if (nrxn == 2) then
reaction_type = "twostep"
else
write(nrxn_string, *) nrxn
reaction_type = trim(nrxn_string) // "-step"
end if
if ( reaction_type == "onestep" ) then
lambda_onestep = pre * exp ( - beta / hrp )
else if ( reaction_type == "twostep" ) then
lambda1_twostep = lambda1
lambda2_twostep = lambda2
beta1_twostep = beta1
hrp_twostep = hrp
else
WRITE(*,*) 'ERROR, UNDEFINED REACTION TYPE ', reaction_type
stop
end if
end subroutine init_chemistry
subroutine update_chemistry (t)
real :: t
real :: factor
real :: relax_duration = 60.
if ( reaction_type == "onestep" ) then
lambda_onestep = pre * exp ( - beta / hrp )
else if ( reaction_type == "twostep" ) then
if (t < relax_duration) then
factor = (relax_duration + t) / relax_duration / 2.
else
factor = 1.
end if
lambda1_twostep = factor * lambda1
lambda2_twostep = factor * lambda2
else
stop
end if
end subroutine update_chemistry
real function rate_1step (yr, theta)
real, intent(in) :: yr
real, intent(in) :: theta
real :: y
real :: t_reduce
y=yr
! if(yr.lt.0.) y=0.
! if(yr.gt.1.) y=1.
t_reduce=theta
! if(theta.lt.0.) t_reduce=0.
! if(theta.gt.1.) t_reduce=1.
if (t_reduce.gt.c_ref) then
rate_1step = pre*y*exp(-ac/(1.+bc*t_reduce))
else if (t_reduce.le.c_cut) then
rate_1step = min_wr
else
rate_1step = &
((refwr-min_wr)*exp(prof_wr*(t_reduce-c_ref)) + min_wr - refwr*exp(prof_wr*(c_cut-c_ref))) &
/ (1.-exp(prof_wr*(c_cut-c_ref)))
endif
end function rate_1step
real function rate1_2step (ya, yx, theta)
real, intent(in) :: ya
real, intent(in) :: yx
real, intent(in) :: theta
real :: y1
real :: y2
real :: t_reduce
y1=ya
if(ya.lt.0.) y1=0.
if(ya.gt.1.) y1=1.
y2=yx
if(yx.lt.0.) y2=0.
if(yx.gt.1.) y2=1.
t_reduce=theta
if(theta.lt.0.) t_reduce=0.
if(theta.gt.1.) t_reduce=1.
rate1_2step = lambda1_twostep * y1 * y2 * &
exp (-(beta1_twostep*(1. - t_reduce))/(1. - hrp_twostep*(1. - t_reduce)))
end function rate1_2step
real function rate2_2step (yx, theta)
real, intent(in) :: yx
real, intent(in) :: theta
real :: y
real :: t_reduce
y=yx
if(yx.lt.0.) y=0.
if(yx.gt.1.) y=1.
t_reduce=theta
if(theta.lt.0.) t_reduce=0.
if(theta.gt.1.) t_reduce=1.
rate2_2step = lambda2_twostep * yx * yx
end function rate2_2step
end module m_chemistry

157
code/m_parameters.f90
View file

@ -1,157 +0,0 @@
module m_parameters
implicit none
! Domain Parameter
integer :: nxp,nyp,nzp
integer :: nx
real :: hxp,hyp,hzp
real :: l_0
real :: hx
! Transport Properties
real :: vis,sc,diff
real :: d_turb
integer :: d_mode
real :: scp,prp,lep,vis0p,rod
real :: lewis, le_a, le_x
! Chemistry Properties
real :: prof_wr,min_wr,min_fsd,min_c,refwr
real :: pre,ac,bc,c_cut,c_ref
real :: lambda1, lambda2, beta1, hrp, beta
character(100) :: reaction_type
! Constants
real, parameter :: pi=3.14159265358979323846d0
! real, parameter :: pi=acos(-1.d0)
real, parameter :: me=1.00e-20
! Flame Control
real :: minf,tar_lo,u0,ctmp,lo_flm=0.
real :: pflame,pflold,oldu
! Solver Control
integer :: nsp, nrxn
integer :: ncyc=0,int_pr
integer :: fctrl_species=1
real :: absolute_tolerence=1e-8
real :: dt,tf,t_now,t_uf,dt_uf
! Input File
LOGICAL :: read_itape, read_stdin
CHARACTER(100) :: itape_name
integer :: table_size
real :: rate_relaxation (2,100)
CHARACTER(LEN=8) :: cdum
INTEGER :: itape=300, otape=301
contains
SUBROUTINE SET_CHEMISTRY
reaction_type = "twostep"
lambda1 = 100000.
lambda2 = lambda1 * 63.6
beta1 = 9.1
hrp = 0.8136765
END SUBROUTINE SET_CHEMISTRY
SUBROUTINE READ_INTRO
IMPLICIT NONE
INTERFACE READ_PARAMETER
SUBROUTINE READ_INT (x)
integer :: x
END SUBROUTINE READ_INT
SUBROUTINE READ_REAL (x)
real :: x
END SUBROUTINE READ_REAL
END INTERFACE
INTEGER :: i
CALL SET_CHEMISTRY
IF (read_stdin) THEN
itape=5
ELSE
OPEN(itape,FILE=itape_name)
END IF
OPEN(otape,FILE='otape')
CALL READ_PARAMETER (nx) ! n grid points
CALL READ_PARAMETER (l_0) ! domain length
CALL READ_PARAMETER (int_pr) ! print interval
CALL READ_PARAMETER (tar_lo) ! target flame location
CALL READ_PARAMETER (dt) ! time step size
CALL READ_PARAMETER (sc) ! Schmidt Number
CALL READ_PARAMETER (vis) ! viscosity
CALL READ_PARAMETER (pre) ! 1-step chemistry
CALL READ_PARAMETER (ac) ! 1-step chemistry
CALL READ_PARAMETER (bc) ! 1-step chemistry
CALL READ_PARAMETER (u0) ! inlet velocity
CALL READ_PARAMETER (tf) ! end time
CALL READ_PARAMETER (dt_uf) ! flame control time step
CALL READ_PARAMETER (ctmp) ! unburned c value
CALL READ_PARAMETER (c_cut) ! reaction rate profile parameter
CALL READ_PARAMETER (c_ref) ! reaction rate profile parameter
CALL READ_PARAMETER (min_wr) ! reaction rate profile parameter
CALL READ_PARAMETER (prof_wr) ! reaction rate profile parameter
CALL READ_PARAMETER (lewis) ! Lewis number
CALL READ_PARAMETER (lambda1) ! 2-step chemistry
CALL READ_PARAMETER (lambda2) ! 2-step chemistry
CALL READ_PARAMETER (beta1) ! 2-step chemistry
CALL READ_PARAMETER (hrp) ! 2-step chemistry
CALL READ_PARAMETER (le_a) ! 2-step chemistry
CALL READ_PARAMETER (le_x) ! 2-step chemistry
CALL READ_PARAMETER (nsp) ! number of species
CALL READ_PARAMETER (nrxn) ! number of reactions
CALL READ_PARAMETER (d_turb) ! turbulent diffusivity
CALL READ_PARAMETER (d_mode) ! diffusivity formula selector
! READ(itape,*) cdum,table_size
! WRITE(otape,*) cdum,table_size
! do i = 1:table_size
! READ(itape,*) cdum,rate_relaxation(1,i),rate_relaxation(2,i)
! WRITE(otape,*) cdum,rate_relaxation(1,i),rate_relaxation(2,i)
! end do
IF (.not.read_stdin) THEN
CLOSE(itape)
END IF
l_0=l_0*pi
hx=l_0/REAL(nx)
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_ref))
END SUBROUTINE READ_INTRO
end module m_parameters

33
code/makefile
View file

@ -1,37 +1,18 @@
BLOCKSIZE?=32 flags = -Wall -O3 -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8 -march=native
flags = -Wall -O3 -cpp -DBLOCKSIZE=$(BLOCKSIZE) -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8 -march=native
ifdef CHECK
flags += -fcheck=all
endif
ifdef DEBUG
flags += -g
endif
compiler = gfortran compiler = gfortran
ex : test.o read_parameter.o ysolve.o m_chemistry.o m_parameters.o m_compact.o ex : test.o ysolve.mod compact.mod
${compiler} -o ex test.o read_parameter.o ysolve.o m_chemistry.o m_parameters.o m_compact.o ${compiler} -o ex test.o ysolve.o Compact.o
test.o : test.f90 ysolve.mod m_compact.mod test.o : test.f90 ysolve.mod
${compiler} -c ${flags} test.f90 ${compiler} -c ${flags} test.f90
read_parameter.o : read_parameter.f90 m_parameters.o ysolve.mod: ysolve.f90 compact.mod
${compiler} -c ${flags} read_parameter.f90
ysolve.o ysolve.mod : ysolve.f90 m_compact.mod m_chemistry.mod m_parameters.mod
${compiler} -c ${flags} ysolve.f90 ${compiler} -c ${flags} ysolve.f90
m_compact.o m_compact.mod : m_compact.f90 compact.mod: Compact.f90
${compiler} -c ${flags} m_compact.f90 ${compiler} -c ${flags} Compact.f90
m_parameters.o m_parameters.mod : m_parameters.f90
${compiler} -c ${flags} m_parameters.f90
m_chemistry.o m_chemistry.mod : m_chemistry.f90 m_parameters.mod
${compiler} -c ${flags} m_chemistry.f90
test: ex test: ex
sh test.sh sh test.sh

32
code/read_parameter.f90
View file

@ -1,32 +0,0 @@
subroutine read_int (x)
use m_parameters
implicit none
integer :: x
READ(itape,*) cdum,x
WRITE(otape,*) cdum,x
end subroutine read_int
subroutine read_real (x)
use m_parameters
implicit none
real :: x
READ(itape,*) cdum,x
WRITE(otape,*) cdum,x
end subroutine read_real
subroutine read_string (x)
use m_parameters
implicit none
character(len=*) :: x
READ(itape,*) cdum,x
WRITE(otape,*) cdum,x
end subroutine read_string

1
code/run.sh
View file

@ -1 +0,0 @@
(TSTAMP=`date -Iseconds`; python sample/itape.py --default sample/4pi-IC1/itape | tee ${TSTAMP}-input.txt | ./ex - | tee ${TSTAMP}-output.txt)

12
code/sample/4pi-IC1/itape
View file

@ -2,24 +2,18 @@
l_0 4. !3[pi] l_0 4. !3[pi]
int_pr 400 int_pr 400
tar_lo 0.60 tar_lo 0.60
dt 0.0005 dt 0.001
sc 0.75 sc 0.75
vis 0.020 !2pi(162)ref !2pi(256)mid !2pi(256)Max !2p(162)Low vis 0.020 !2pi(162)ref !2pi(256)mid !2pi(256)Max !2p(162)Low
pre 2.10E+4 !185.75 !185.75 !316.88 !440.5 !102.69 pre 2.10E+4 !185.75 !185.75 !316.88 !440.5 !102.69
ac 26.7 !22.250 !22.25 !22.56 !22.3 !22.39 ac 26.7 !22.250 !22.25 !22.56 !22.3 !22.39
bc 3. !7.000 !7.0 !7.0 !7.0 !7.0 bc 3. !7.000 !7.0 !7.0 !7.0 !7.0
u0 0.21 u0 0.21
tf 200. tf 80.
dt_uf 0.4 dt_uf 0.4
ctmp 0. !1.0e-14 ctmp 0. !1.0e-14
c_cut 0.001 !0.012 ! c < c_cut -> wrate = 0. c_cut 0.001 !0.012 ! c < c_cut -> wrate = 0.
c_ref 0.01 !0.003 ! c_ref 0.01 !0.003 !
min_wr 0. ! 5.0e-14 min_wr 0. ! 5.0e-14
prof_wr 1.0 prof_wr 1.0
lewis 1.0
lambda1 6.25E+4 7.7E+4 5E+4
lambda2 318E+4 445.2E+4 489.72E+4
beta1 5
hrp 0.75
le_a 1.
le_x 0.20

89
code/sample/itape.py
View file

@ -1,89 +0,0 @@
#!/bin/python
import os, sys
import argparse
from pprint import pprint
parser = argparse.ArgumentParser()
file_parser = argparse.ArgumentParser()
# Default values for all parameters
raw_defaults = (
'''\
nx 512
l_0 4.
int_pr 400
tar_lo 0.60
dt 0.0005
sc 0.75
vis 0.020
pre 2.10E+4
ac 26.7
bc 3.
u0 0.21
tf 100.
dt_uf 0.4
ctmp 0.
c_cut 0.001
c_ref 0.01
min_wr 0.
prof_wr 1.0
lewis 1.0
lambda1 6.25E+4
lambda2 318E+4
beta1 5
hrp 0.75
le_a 1.
le_x 0.20
nsp 2
nrxn 1
d_turb 0
d_mode 0
'''
).split()
ordered_args = raw_defaults[::2]
ordered_values = []
for v in raw_defaults[1::2]:
try:
ordered_values.append(int(v))
except ValueError:
ordered_values.append(float(v))
param_dict = dict(zip(ordered_args, ordered_values))
file_parser.add_argument("--default", help="file containing override default values")
file_args, others = file_parser.parse_known_args()
if file_args.default is not None:
with open(file_args.default) as df:
try:
for line in df:
k, v = line.split()[:2]
if k in param_dict:
param_dict[k] = v
except ValueError:
pass
for k, v in param_dict.iteritems():
parser.add_argument("--"+k, default=v)
args = parser.parse_args(others)
params = vars(args)
line_new = '{:>12} {}'
for a in ordered_args:
print line_new.format(a, params[a])
# for i in range(15,60,5):
# with open("itape", 'w') as itape:
# itape.write(format.format(float(i)/10.0))
# os.system("../../ex &>"+str(i)+".txt")

7
code/test.f90
View file

@ -6,12 +6,5 @@
PROGRAM test PROGRAM test
USE ysolve USE ysolve
IMPLICIT NONE IMPLICIT NONE
CALL parse
CALL init_solver
CALL solve CALL solve
CALL finalize_solver
END PROGRAM test END PROGRAM test

4
code/test.sh
View file

@ -1,4 +1,4 @@
cd sample/4pi-IC1 cd sample/4pi-IC1
ls ls
../../ex itape ../../ex
diff -y sfield.dat sfield.diff diff sfield.dat sfield.diff

785
code/ysolve.f90
View file

@ -1,202 +1,76 @@
MODULE ysolve MODULE ysolve
USE m_compact USE Compact
USE m_parameters
USE m_chemistry
IMPLICIT NONE IMPLICIT NONE
PRIVATE
REAL, DIMENSION(:), ALLOCATABLE :: u REAL, PARAMETER :: pi=3.14159265358979323846
REAL, DIMENSION(:), ALLOCATABLE :: inletbc REAL :: hx,dt,vis,sc,diff,pre,ac,bc,tf,t_now,t_uf,dt_uf
REAL, DIMENSION(:,:), ALLOCATABLE :: y1,y2,yf,yold 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 REAL, DIMENSION(:), ALLOCATABLE :: uxt,duxt
REAL, DIMENSION(:), ALLOCATABLE :: dm
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 CONTAINS
!------------------------------------------------------------------------ !------------------------------------------------------------------------
SUBROUTINE parse
CHARACTER(100) :: num1char
! Good. One command line argument provided.
IF(COMMAND_ARGUMENT_COUNT().EQ.1)THEN
CALL GET_COMMAND_ARGUMENT(1,num1char)
! Bad. No command line argument provided. Print help message.
ELSE IF(COMMAND_ARGUMENT_COUNT().EQ.0)THEN
WRITE(*,*)'ERROR, NO COMMAND-LINE ARGUMENT'
num1char = "-h"
! Bad. Multiple command line arguments provided. Print help message.
ELSE
WRITE(*,*)'ERROR, TOO MANY COMMAND-LINE ARGUMENTS( > 2 )'
num1char = "-h"
ENDIF
! Option or STDIN
IF(num1char(1:1) == "-") THEN
if(num1char=="-h") then
write(*,'(a)') "usage: ex [-h] -|input_file"
write(*,'(a)') ""
write(*,'(a)') "positional arguments:"
write(*,'(a)') " - read from std_in"
write(*,'(a)') " input_file input file name"
write(*,'(a)') ""
write(*,'(a)') "optional arguments:"
write(*,'(a)') " -h show this help message and exit"
stop
else if (num1char=="-") then
read_stdin = .true.
else
WRITE(*,*)'ERROR, UNSUPPORTED OPTION ', trim(num1char), '. STOPPING'
STOP
end if
END IF
IF(read_stdin)THEN
WRITE(*,*) "Read from STDIN"
ELSE
itape_name = num1char
WRITE(*,*) "Read from file " // num1char
END IF
END SUBROUTINE parse
!------------------------------------------------------------------------
REAL FUNCTION residual (x0, x1)
REAL, DIMENSION(:,:) :: x0, x1
residual = sum(abs(x0-x1))
END FUNCTION residual
LOGICAL FUNCTION converged (x0, x1)
REAL, DIMENSION(:,:) :: x0, x1
REAL :: r = 0.
r = residual(x0 ,x1)
converged = r < absolute_tolerence
END FUNCTION converged
!------------------------------------------------------------------------
SUBROUTINE solve SUBROUTINE solve
IF (d_mode == 0) THEN INTEGER :: i,j,k,savenum
CALL solve_with_diffusivity (update_dm)
ELSE IF (d_mode == 1) THEN
CALL solve_with_diffusivity (update_dm_cd)
ELSE IF (d_mode == 2) THEN
CALL solve_with_diffusivity (update_dm_sutherland)
ELSE
WRITE(*,*)'ERROR, UNSUPPORTED DIFFUSIVITY ', d_mode, '. STOPPING'
STOP
END IF
END SUBROUTINE solve
SUBROUTINE solve_with_diffusivity (calc_diff)
INTERFACE
SUBROUTINE calc_diff(r1)
REAL, INTENT(IN), DIMENSION(:,:) :: r1
END SUBROUTINE calc_diff
END INTERFACE
IF (reaction_type == "onestep") THEN
CALL solve_ (fonestep, SET_IC_ONESTEP, calc_diff)
ELSE IF (reaction_type == "twostep") THEN
CALL solve_ (fns, SET_IC_TWOSTEP, calc_diff)
ELSE
WRITE(*,*)'ERROR, UNSUPPORTED REACTION ', trim(reaction_type), '. STOPPING'
STOP
END IF
END SUBROUTINE solve_with_diffusivity
SUBROUTINE solve_ (rhs, SET_IC, calc_diff)
INTERFACE
SUBROUTINE rhs(r1,f)
REAL, INTENT(IN), DIMENSION(:,:) :: r1
REAL, INTENT(OUT), DIMENSION(:,:) :: f
END SUBROUTINE rhs
SUBROUTINE SET_IC
END SUBROUTINE SET_IC
SUBROUTINE calc_diff(r1)
REAL, INTENT(IN), DIMENSION(:,:) :: r1
END SUBROUTINE calc_diff
END INTERFACE
INTEGER :: i
REAL :: pflame,pflold,delf=0. REAL :: pflame,pflold,delf=0.
REAL :: residue = 0.
CALL READ_INTRO
CALL ludcmp(nx,5,5,1,0,0)
CALL SET_IC CALL SET_IC
pflame=0. t_now=0.; t_uf=0.
DO i=1,nx
pflame=pflame+y1(i,fctrl_species)*hx
ENDDO
t_now=0.
t_uf=0.
DO 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
IF (t_now.ge.tf) EXIT ! 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(*,*)
residue = residual(yold, y1) ENDIF
IF (converged(yold, y1)) EXIT
ncyc=ncyc+1 ncyc=ncyc+1
t_uf=t_uf+dt t_uf=t_uf+dt
t_now=t_now+dt t_now=t_now+dt
yold = y1
CALL update_chemistry(t_now)
CALL SET_BC
CALL RK4(rhs, calc_diff)
IF (t_uf.ge.dt_uf) THEN
pflold=pflame
pflame=0.
DO i=1,nx
pflame=pflame+y1(i,fctrl_species)*hx
uxt(i)= 1. - y1(i,fctrl_species)
ENDDO
CALL dfnonp(nx,hx,uxt,duxt,1,1)
delf=1./MAXVAL(ABS(duxt))
oldu=u(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,a6,e10.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),' // L_f:',delf,' // R:',residue
ENDIF
IF (MOD(ncyc,int_pr).eq.0) THEN IF (MOD(ncyc,int_pr).eq.0) THEN
! WRITE(*,'(a2,f8.3,a9,f10.7,a11,i6,a7,f9.5)') & ! WRITE(*,'(a2,f8.3,a9,f10.7,a11,i6,a7,f9.5)') &
! 'T:',t_now,' // dT:',dt,' // NCYC:',ncyc,' // U:',u(1,1,1) ! 'T:',t_now,' // dT:',dt,' // NCYC:',ncyc,' // U:',u(1,1,1)
ENDIF ENDIF
ENDDO ENDDO
! CALL write_sd CALL write_sd
CALL write_pre CALL write_pre
CALL save_final_field CALL save_final_field
@ -204,12 +78,12 @@
WRITE(*,*) 'Fin.' WRITE(*,*) 'Fin.'
WRITE(*,*) WRITE(*,*)
END SUBROUTINE solve_ END SUBROUTINE solve
SUBROUTINE write_sd SUBROUTINE write_sd
REAL :: c,yr,theta,wrate,dely,sdr,sdd,sd,uu,onelw,dd REAL :: c,yr,wrate,dely,sdr,sdd,sd,uu,onelw,dd
INTEGER :: i,j,k,nd INTEGER :: i,j,k,nd
REAL, DIMENSION(2,nx) :: ux,dux,d2ux REAL, DIMENSION(1,nx) :: ux,dux,d2ux
REAL, DIMENSION(10,nx) :: sav REAL, DIMENSION(10,nx) :: sav
sav=0. sav=0.
@ -223,22 +97,21 @@
WRITE(500,*) ' "Sdd","Sd","(1/C)/(dC/dx)","DIV(rho*Dmu*Gra(C))"' WRITE(500,*) ' "Sdd","Sd","(1/C)/(dC/dx)","DIV(rho*Dmu*Gra(C))"'
DO i=1,nx DO i=1,nx
ux(1,i)=y1(i,1) ! Yr ux(1,i)=y1(i,1,1) ! Yr
ux(2,i)=y1(i,2) ! T
IF (ux(1,i).gt.1.) ux(1,i)=1. IF (ux(1,i).gt.1.) ux(1,i)=1.
ENDDO ENDDO
nd=2 nd=1
CALL dfnonp(nx,hx,ux(:,:),dux(:,:),nd,1) CALL dfnonp(nx,hx,ux(1,:),dux(1,:),nd,1)
CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),nd,1) nd=1
CALL d2fnonp(nx,hx,ux(1,:),d2ux(1,:),nd,1)
DO i=1,nx DO i=1,nx
yr=ux(1,i) yr=ux(1,i)
theta=ux(2,i)
c=1.-yr c=1.-yr
IF (c.lt.0.) c=0. IF (c.lt.0.) c=0.
wrate=pre*yr*exp(-ac/(1.+bc*(theta))) !wrate wrate=pre*yr*exp(-ac/(1.+bc*(1.-yr))) !wrate
! IF (c.le.c_cut) THEN ! IF (c.le.c_cut) THEN
! wrate=min_wr ! wrate=min_wr
@ -246,22 +119,22 @@
! ((exp((c-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr ! ((exp((c-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr
! ENDIF ! ENDIF
IF (theta.le.c_ref) THEN IF (c.le.c_ref) THEN
wrate=min_wr wrate=min_wr
IF (theta.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(theta-c_ref))+ & 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))) min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
ENDIF ENDIF
dely=ABS(dux(1,i)) dely=ABS(dux(1,i))
sdr=wrate/dely sdr=wrate/dely
sdd=-dm(i)*d2ux(1,i)/dely sdd=-diff*d2ux(1,i)/dely
IF (dely.eq.0.) THEN IF (dely.eq.0.) THEN
sdr=0.; sdd=0. sdr=0.; sdd=0.
ENDIF ENDIF
sd=sdr+sdd sd=sdr+sdd
uu=u(1) uu=u(1,1,1)
onelw=(-dux(1,i))/c onelw=(-dux(1,i))/c
dd=-dm(i)*d2ux(1,i) dd=-diff*d2ux(1,i)
if (c.eq.0.) onelw=0. if (c.eq.0.) onelw=0.
sav(1,i)=sav(1,i)+yr sav(1,i)=sav(1,i)+yr
@ -282,63 +155,48 @@
END SUBROUTINE write_sd END SUBROUTINE write_sd
SUBROUTINE save_final_field SUBROUTINE save_final_field
INTEGER :: i INTEGER :: i,j,k
OPEN (305,FILE='sfield.bin',form='unformatted',status='unknown') OPEN (305,FILE='sfield.bin',form='unformatted',status='unknown')
DO i=1,nx DO i=1,nx
WRITE (305) y1(i,1) WRITE (305) y1(i,1,1)
ENDDO ENDDO
CLOSE (305) CLOSE (305)
OPEN (305,FILE='sfield.dat') OPEN (305,FILE='sfield.dat')
DO i=1,nx DO i=1,nx
WRITE (305,*) (/ i*hx, y1(i,:) /) WRITE (305,'(e30.20)') y1(i,1,1)
ENDDO ENDDO
CLOSE (305) CLOSE (305)
END SUBROUTINE save_final_field END SUBROUTINE save_final_field
SUBROUTINE write_pre SUBROUTINE write_pre
REAL :: theta,yr,c,dy,maxdy=0.,del_f REAL :: yr,c,dy,maxdy=0.,del_f
REAL :: S_L=0.,wrate,wrate1,wrate2 REAL :: S_L=0.,wrate
REAL :: ya,yx
INTEGER :: i INTEGER :: i
REAL, DIMENSION(1,nx) :: ux, dux REAL, DIMENSION(1,nx) :: ux, dux
if ( reaction_type == "onestep" ) then
DO i=1,nx DO i=1,nx
yr=y1(i,1) yr=y1(i,1,1)
theta=y1(i,2) c=1.-yr
ux(1,i)=yr
ux(1,i)=theta wrate=pre*yr*exp(-ac/(1.+bc*(1.-yr))) !wrate
! IF (c.le.c_cut) THEN
wrate=rate_1step(yr, theta) ! 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 S_L=S_L+wrate*hx
ENDDO ENDDO
else if ( reaction_type == "twostep" ) then
DO i=1,nx
ya=y1(i,1)
yx=y1(i,2)
theta=y1(i,3)
ux(1,i)=theta
wrate1=rate1_2step(ya, yx, theta)
wrate2=rate2_2step(yx, theta)
S_L=S_L+wrate1*hx
ENDDO
else
WRITE(*,*) 'ERROR, UNDEFINED REACTION TYPE ', reaction_type
stop
end if
WRITE(*,'(a31,e14.8)') ' INTEGRAL( Wrate x dx ) => Sc :',S_L WRITE(*,'(a31,e14.8)') ' INTEGRAL( Wrate x dx ) => Sc :',S_L
CALL dfnonp(nx,hx,ux(1,:),dux(1,:),1,1) CALL dfnonp(nx,hx,ux(1,:),dux(1,:),1,1)
@ -349,337 +207,252 @@
ENDDO ENDDO
del_f=1./maxdy del_f=1./maxdy
WRITE(*,'(a13,e14.8,a25,e14.8)') ' Grid size : ',hx,' / Laminar flame speed : ',u(1) 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)') & WRITE(*,'(a19,e14.8,a3,f9.5,a20)') &
' Flame thickness : ',del_f,' / ',del_f/hx,' grids in the flame.' ' Flame thickness : ',del_f,' / ',del_f/hx,' grids in the flame.'
WRITE(*,*) WRITE(*,*)
END SUBROUTINE write_pre END SUBROUTINE write_pre
SUBROUTINE SET_BC SUBROUTINE SET_BC
INTEGER :: i
DO i = 1, nsp y1(1,1,1)=1.-ctmp
y1(1,i) = inletbc(i)
END DO
END SUBROUTINE SET_BC END SUBROUTINE SET_BC
SUBROUTINE BASE_FLAME_PROFILE(x) SUBROUTINE SET_IC
INTEGER :: i, ifl, si INTEGER :: i, ifl, si
REAL :: x(nx) REAL :: xi
! initial flame thickness (0.2 pi n grids)
ifl=INT(2.0*pi*0.1/hx) ifl=INT(2.0*pi*0.1/hx)
! initialize velocity field u=u0; si=INT(nx*(1.-tar_lo))
u=u0;
! initial flame center (grid index)
si=INT(nx*(1.-tar_lo))
! initialize Yr field
DO i=1,nx DO i=1,nx
IF(i< nx-(si+ifl/2)) THEN IF(i< nx-(si+ifl/2)) THEN
x(i)=0.+ctmp xi=0.+ctmp
ELSE IF(i> nx-(si-ifl/2)) THEN ELSE IF(i> nx-(si-ifl/2)) THEN
x(i)=1. xi=1.
ELSE ELSE
x(i)=0.5+REAL(i-nx+si)/REAL(ifl) xi=0.5+REAL(i-nx+si)/REAL(ifl)
ENDIF ENDIF
END DO
! initialize species field ! y1(i,1,1)=(1.-xi)*1.
! max_ysum = maxval(y1(:,1) + y1(:,2)) y1(i,1,1)=(1.-xi) ! reactant mass fraction
! y1(:,1) = y1(:,1) / maxval(y1(:,1) + y1(:,2)) ENDDO
! y1(:,2) = y1(:,2) / maxval(y1(:,1) + y1(;,2))
! y1(:,1) = y1(:,1) / (y1(1,1) + y1(1,2))
! y1(:,2) = y1(:,2) / (y1(1,1) + y1(1,2))
END SUBROUTINE BASE_FLAME_PROFILE pflame=0.
DO i=1,nx
pflame=pflame+y1(i,1,1)*hx
ENDDO
SUBROUTINE SET_IC_ONESTEP END SUBROUTINE SET_IC
INTEGER :: i
REAL :: xi
REAL :: x(nx)
CALL BASE_FLAME_PROFILE(x) SUBROUTINE READ_INTRO
CHARACTER(LEN=8) :: cdum
INTEGER :: itape=300,otape=301,ierr
DO i=1,nx
xi = x(i)
y1(i,1)=(1.-xi) ! reactant mass fraction OPEN(itape,FILE='itape')
y1(i,2)=(xi) ! reactant mass fraction OPEN(otape,FILE='otape')
ENDDO READ(itape,*) cdum,nx
END SUBROUTINE SET_IC_ONESTEP 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
SUBROUTINE SET_IC_TWOSTEP CLOSE(itape)
INTEGER :: i
REAL :: xi
REAL :: x(nx)
CALL BASE_FLAME_PROFILE(x) ! 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
DO i=1,nx ! refwr=pre*1.*exp(-ac/(1.+bc*c_cut))
xi = x(i) ! minf=exp((c_ref-c_cut)*prof_wr)
y1(i,1)=(1.-xi) ! reactant mass fraction refwr=pre*1.*exp(-ac/(1.+bc*c_ref))
l_0=l_0*pi
y1(i,2) = (1./2.) * (lambda1/lambda2) * y1(i,1) * & ALLOCATE(u(nx,1,1),STAT=ierr) ; u=0.
exp (-(beta1*(1. - xi))/(1. - hrp*(1. - xi))) ALLOCATE(y1(nx,1,1),STAT=ierr) ; y1=0.
ALLOCATE(y2(nx,1,1),STAT=ierr) ; y2=0.
y1(i,3) = xi ALLOCATE(yf(nx,1,1),STAT=ierr) ; yf=0.
ENDDO
END SUBROUTINE SET_IC_TWOSTEP
SUBROUTINE init_solver
INTEGER :: ierr
CALL READ_INTRO
CALL init_chemistry
ALLOCATE( u(nx),STAT=ierr) ; u=0.
ALLOCATE( inletbc(nsp),STAT=ierr) ; inletbc=0.
ALLOCATE( y1(nx,nsp),STAT=ierr) ; y1=0.
ALLOCATE( y2(nx,nsp),STAT=ierr) ; y2=0.
ALLOCATE( yf(nx,nsp),STAT=ierr) ; yf=0.
ALLOCATE(yold(nx,nsp),STAT=ierr) ; yold=0.
ALLOCATE(uxt(nx),STAT=ierr) ; uxt=0. ALLOCATE(uxt(nx),STAT=ierr) ; uxt=0.
ALLOCATE(duxt(nx),STAT=ierr) ; duxt=0. ALLOCATE(duxt(nx),STAT=ierr) ; duxt=0.
ALLOCATE(dm(nx),STAT=ierr) ; dm=diff
if ( reaction_type == "onestep" ) then END SUBROUTINE READ_INTRO
inletbc(1) = 1.
else if ( reaction_type == "twostep" ) then
inletbc(1) = 1.
inletbc(2) = 0.
inletbc(3) = 0.
else
WRITE(*,*) 'ERROR, UNDEFINED REACTION TYPE ', reaction_type
stop
end if
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
CALL ludcmp(nx,100,100,1,0,0)
END SUBROUTINE init_solver
SUBROUTINE finalize_solver
DEALLOCATE( u)
DEALLOCATE(inletbc)
DEALLOCATE(y1)
DEALLOCATE(y2)
DEALLOCATE(yf)
DEALLOCATE(yold)
DEALLOCATE(uxt)
DEALLOCATE(duxt)
DEALLOCATE(dm)
END SUBROUTINE finalize_solver
SUBROUTINE RK4(rhs, calc_diff)
INTERFACE
SUBROUTINE rhs(r1,f)
REAL, INTENT(IN), DIMENSION(:,:) :: r1
REAL, INTENT(OUT), DIMENSION(:,:) :: f
END SUBROUTINE rhs
SUBROUTINE calc_diff(r1)
REAL, INTENT(IN), DIMENSION(:,:) :: r1
END SUBROUTINE calc_diff
END INTERFACE
CALL substep(1,y1,y1,y2,yf)
CALL substep(2,y1,y2,y1,yf)
CALL substep(3,y2,y1,y2,yf)
CALL substep(4,y1,y2,y1,yf)
CALL substep(5,y2,y1,y2,yf)
CONTAINS
SUBROUTINE substep(istage, ri,r1,r2,f)
REAL, DIMENSION(5), PARAMETER :: a=(/ 970286171893.d0/4311952581923., &
6584761158862.d0/12103376702013., &
2251764453980.d0/15575788980749., &
26877169314380.d0/34165994151039., &
0.d0 /), &
b=(/ 1153189308089.d0/22510343858157., &
1772645290293.d0/4653164025191., &
-1672844663538.d0/4480602732383., &
2114624349019.d0/3568978502595., &
5198255086312.d0/14908931495163. /)
INTEGER :: istage
REAL, INTENT(INOUT),DIMENSION(:,:) :: ri,r1,r2
REAL, INTENT(OUT),DIMENSION(:,:) :: f
REAL :: at,bt REAL :: at,bt
CALL calc_diff(ri) CALL fns(ri,f)
CALL rhs(ri,f)
IF(istage<5) THEN IF(istage<5) THEN
at=a(istage)*dt at=a(istage)*dt
bt=(b(istage)-a(istage))*dt bt=(b(istage)-a(istage))*dt
r1=r1+at*f DO k=1,1 ! nz
r2=r1+bt*f 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 ELSE
bt=b(istage)*dt bt=b(istage)*dt
r1=r1+bt*f 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 ENDIF
END SUBROUTINE substep END SUBROUTINE substep
!------------------------------------------------------------------------
END SUBROUTINE RK4
SUBROUTINE update_dm(r1)
REAL, INTENT(IN),DIMENSION(:,:) :: r1
INTEGER :: i
DO i=1,nx
dm(i) = diff + d_turb
ENDDO
END SUBROUTINE update_dm
SUBROUTINE update_dm_cd(r1)
REAL, INTENT(IN),DIMENSION(:,:) :: r1
REAL :: conv_rxn_boundary
REAL :: beta, tutb
INTEGER :: i
tutb = 1. / (bc+1.)
! \beta = Ta/Tb * (Tb-Tu)/Tb
beta = ac * tutb * (1. - tutb)
conv_rxn_boundary = 1./beta
DO i=1,nx
IF (r1(i,fctrl_species) > conv_rxn_boundary) THEN
dm(i) = diff + d_turb
ELSE
dm(i) = diff
END IF
ENDDO
END SUBROUTINE update_dm_cd
SUBROUTINE update_dm_sutherland(r1)
REAL, INTENT(IN),DIMENSION(:,:) :: r1
INTEGER :: i
DO i=1,nx
dm(i) = diff * diffusivity_sutherland(r1(i,nsp)) + d_turb
ENDDO
END SUBROUTINE update_dm_sutherland
SUBROUTINE fonestep(r1,f)
REAL, INTENT(IN),DIMENSION(:,:) :: r1
REAL, INTENT(OUT),DIMENSION(:,:) :: f
REAL, DIMENSION(nsp,nx) :: ux, dux, d2ux
REAL, DIMENSION(nx) :: dxdm
INTEGER :: i
REAL :: wrate,Ly,Dy,Lt,Dt
! x-direction
DO i=1,nx
ux(1,i)=r1(i,1) ! Y
ux(2,i)=r1(i,2) ! T
ENDDO
CALL dfnonp(nx,hx,ux(:,:),dux(:,:),nsp,1)
CALL dfnonp(nx,hx,dm,dxdm,1,1)
CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),nsp,1)
DO i=1,nx
wrate=rate_1step(ux(1,i), ux(2,i))
! - u*dY/dx + D*d2Y/d2x
f(i,1) = - ( u(i)*dux(1,i) ) + dm(i) * d2ux(1,i) + dxdm(i) * dux(1,i) - wrate
! - u*dY/dx + D*d2Y/d2x
f(i,2) = - ( u(i)*dux(2,i) ) + dm(i) * d2ux(2,i) + dxdm(i) * dux(2,i) + wrate
! Boundary conditions
IF (i.eq.nx) THEN
f(nx,1) = -wrate - u(nx)*dux(1,nx)
f(nx,2) = wrate - u(nx)*dux(2,nx)
ENDIF
ENDDO
! Boundary conditions
f(1,1)=0.
f(1,2)=0.
END SUBROUTINE fonestep
SUBROUTINE fns(r1,f) SUBROUTINE fns(r1,f)
REAL, INTENT(IN),DIMENSION(:,:) :: r1 REAL, INTENT(IN),DIMENSION(:,:,:) :: r1
REAL, INTENT(OUT),DIMENSION(:,:) :: f REAL, INTENT(OUT),DIMENSION(:,:,:) :: f
REAL, DIMENSION(3,nx) :: ux, dux, d2ux REAL, DIMENSION(3,nx) :: ux,dux,d2ux
INTEGER :: i INTEGER :: i,j,k
REAL :: Ly,Dy,Lt,Dt REAL :: wrate,Ly,Dy
REAL :: wrate1, wrate2
DO k=1,1 !nz
! x-direction ! x-direction
DO i=1,nx DO j=1,1 !ny
ux(1,i)=r1(i,1) ! Y DO i=1,nx
ux(2,i)=r1(i,2) ! Y ux(1,i)=u(i,j,k)*r1(i,j,k) ! u*Y
ux(3,i)=r1(i,3) ! T 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 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
CALL dfnonp(nx,hx,ux(:,:),dux(:,:),3,1) ! 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
CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),3,1) ! Boundary conditionS
f(1,1,1)=0.
DO i=1,nx Dy=Ly
wrate1=rate1_2step(ux(1,i), ux(2,i), ux(3,i)) f(nx,1,1)=f(nx,1,1)-Dy
wrate2=rate2_2step(ux(2,i), ux(3,i))
! - u*dY/dx + D*d2Y/d2x
f(i,1) = - ( u(i)*dux(1,i) ) + (diff/le_a) * d2ux(1,i) - wrate1
! - u*dY/dx + D*d2Y/d2x
f(i,2) = - ( u(i)*dux(2,i) ) + (diff/le_x) * d2ux(2,i) + wrate1 - 2.*wrate2
! - u*dT/dx + D*d2T/d2x
f(i,3) = - ( u(i)*dux(3,i) ) + (diff) * d2ux(3,i) + 2.*wrate2
! Boundary conditions
IF (i.eq.nx) THEN
f(nx,1) = -wrate1 - u(nx)*dux(1,nx)
f(nx,2) = wrate1 - 2.*wrate2 - u(nx)*dux(2,nx)
f(nx,3) = 2.*wrate2 - u(nx)*dux(3,nx)
ENDIF
ENDDO
! Boundary conditions
f(1,1)=0.
f(1,2)=0.
f(1,3)=0.
END SUBROUTINE fns END SUBROUTINE fns
!------------------------------------------------------------------------
REAL FUNCTION diffusivity_sutherland(c)
REAL, INTENT(IN) :: c
REAL :: theta, As, Ts, T0, T1
REAL, PARAMETER :: rvis=5.0d0
T0 = 1.0
T1 = (1.0 + bc)
Ts = (rvis*T1 - (T1**(3./2.))) / (T1**(3./2.) - rvis)
As = (T0 + Ts)
theta = (1.0 + bc * c)
diffusivity_sutherland = As * sqrt(theta) / (1. + Ts/theta)
END FUNCTION diffusivity_sutherland
END MODULE ysolve END MODULE ysolve