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:
- cd ${CI_PROJECT_DIR}
- cd code
- make clean && make
- cd sample/4pi-IC1
- ../../ex
- make test

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

@ -1,13 +1,14 @@
MODULE m_compact
MODULE Compact
IMPLICIT NONE
REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, &
lyf,lys,wyf,wys, &
lzf,lzs,wzf,wzs
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(KIND=8), PARAMETER :: ezero = 1.0e-14
REAL, PARAMETER :: ezero = 1.0e-14
PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp
CONTAINS
@ -19,22 +20,6 @@
nxc=nx
nyc=ny
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'
! xp, yp, zp = 0 : periodic
@ -47,6 +32,9 @@
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)
@ -58,6 +46,9 @@
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)
@ -69,140 +60,19 @@
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'
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)
ELSE
call nonp_lud(3,nzc)
ENDIF
END SUBROUTINE ludcmp_calculate
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
! CALL x_lud
! CALL yz_lud
END SUBROUTINE ludcmp
SUBROUTINE nonp_lud(xyz,xx)
INTEGER :: i,xyz,xx
REAL(KIND=8), DIMENSION(xx) :: aa
REAL, DIMENSION(xx) :: aa
aa=3.
aa(1)=0.5 ; aa(2)=4.
aa(xx-1)=4. ; aa(xx)=0.5
@ -222,7 +92,7 @@
SUBROUTINE p_lud(xyz,xx)
INTEGER :: i,xyz,xx
REAL(KIND=8) :: a
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
@ -235,23 +105,21 @@
SUBROUTINE stdlu(a,n,l)
INTEGER :: n
REAL(KIND=8), INTENT(IN) :: a(n)
REAL(KIND=8), INTENT(OUT) :: l(n)
REAL(KIND=8) :: d
REAL :: a(n),l(n)
REAL :: d
INTEGER :: i
l(1)=1.0d0/a(1)
l(1)=1.0/a(1)
DO i=2,n
d=a(i)-l(i-1)
l(i)=1.0d0/d
l(i)=1.0/d
ENDDO
END SUBROUTINE stdlu
SUBROUTINE ptdlu(a,n,l,w)
INTEGER :: n
REAL(KIND=8), INTENT(OUT) :: a
REAL(KIND=8), INTENT(OUT) :: l(n),w(n)
REAL :: a,l(n),w(n)
INTEGER :: i
REAL(KIND=8) :: aa(n),d
REAL :: aa(n),d
DO i=1,n-1
aa(i)=a
@ -273,19 +141,18 @@
l(n)=1./d
END SUBROUTINE ptdlu
SUBROUTINE rhs1np(n,h,x,dx,nd)
INTEGER,INTENT(IN) :: n,nd
REAL(KIND=8),INTENT(IN) :: h
REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
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(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
r2=1.d0/12.d0
r1=7./3.
r2=1./12.
r3=3.
a=-1.25
b=1.
@ -311,20 +178,6 @@
dx(j,i)=h1*(r1*t1+r2*t2)
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.2) CALL tdslv(dx,n,lyf,nd) ! y-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)
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
REAL,INTENT(IN) :: h
REAL,INTENT(IN),DIMENSION(nd,n) :: x
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j
REAL(KIND=8) :: r1,r2,h1
! print *, "dfnonp received (nd,n)", nd, n
REAL :: r1,r2,h1
h1=1./h
r1=7./3.
@ -371,10 +221,10 @@
SUBROUTINE ptdslv(r,n,l,w,nd)
INTEGER,INTENT(IN) :: n,nd
REAL(KIND=8),INTENT(INOUT),DIMENSION(nd,n) :: r
REAL(KIND=8),INTENT(IN),DIMENSION(n) :: l,w
REAL,INTENT(INOUT),DIMENSION(nd,n) :: r
REAL,INTENT(IN),DIMENSION(:) :: l,w
INTEGER i,j
REAL(KIND=8), DIMENSION(nd) :: sum
REAL, DIMENSION(nd) :: sum
DO j=1,nd
sum(j)=w(1)*r(j,1)
r(j,1)=r(j,1)*l(1)
@ -399,13 +249,11 @@
SUBROUTINE d2fp(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
REAL,INTENT(IN) :: h
REAL,INTENT(IN),DIMENSION(nd,n) :: x
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j
REAL(KIND=8) :: h2,r1,r2,t1,t2
REAL :: h2,r1,r2,t1,t2
h2=1./(h*h)
r1=6.
r2=3./8.
@ -464,10 +312,10 @@
SUBROUTINE tdslv(r,n,l,nd)
INTEGER,INTENT(IN) :: n,nd
REAL(KIND=8),INTENT(INOUT),DIMENSION(nd,n) :: r
REAL(KIND=8),INTENT(IN),DIMENSION(n) :: l
REAL,INTENT(INOUT),DIMENSION(nd,n) :: r
REAL,INTENT(IN),DIMENSION(:) :: l
INTEGER i,j
REAL(KIND=8) t1
REAL t1
DO j=1,nd
r(j,1)=r(j,1)*l(1)
ENDDO
@ -486,12 +334,11 @@
SUBROUTINE d2fnonp(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
REAL,INTENT(IN) :: h
REAL,INTENT(IN),DIMENSION(nd,n) :: x
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
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)
r1=6.
@ -536,4 +383,4 @@
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 -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
flags = -Wall -O3 -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8 -march=native
compiler = gfortran
ex : test.o read_parameter.o ysolve.o m_chemistry.o m_parameters.o m_compact.o
${compiler} -o 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 ysolve.o Compact.o
test.o : test.f90 ysolve.mod m_compact.mod
test.o : test.f90 ysolve.mod
${compiler} -c ${flags} test.f90
read_parameter.o : read_parameter.f90 m_parameters.o
${compiler} -c ${flags} read_parameter.f90
ysolve.o ysolve.mod : ysolve.f90 m_compact.mod m_chemistry.mod m_parameters.mod
ysolve.mod: ysolve.f90 compact.mod
${compiler} -c ${flags} ysolve.f90
m_compact.o m_compact.mod : m_compact.f90
${compiler} -c ${flags} m_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
compact.mod: Compact.f90
${compiler} -c ${flags} Compact.f90
test: ex
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]
int_pr 400
tar_lo 0.60
dt 0.0005
dt 0.001
sc 0.75
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
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
u0 0.21
tf 200.
tf 80.
dt_uf 0.4
ctmp 0. !1.0e-14
c_cut 0.001 !0.012 ! c < c_cut -> wrate = 0.
c_ref 0.01 !0.003 !
min_wr 0. ! 5.0e-14
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
USE ysolve
IMPLICIT NONE
CALL parse
CALL init_solver
CALL solve
CALL finalize_solver
END PROGRAM test

4
code/test.sh
View file

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

803
code/ysolve.f90
View file

@ -1,202 +1,76 @@
MODULE ysolve
USE m_compact
USE m_parameters
USE m_chemistry
USE Compact
IMPLICIT NONE
REAL, DIMENSION(:), ALLOCATABLE :: u
REAL, DIMENSION(:), ALLOCATABLE :: inletbc
REAL, DIMENSION(:,:), ALLOCATABLE :: y1,y2,yf,yold
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
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
!------------------------------------------------------------------------
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
IF (d_mode == 0) THEN
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
INTEGER :: i,j,k,savenum
REAL :: pflame,pflold,delf=0.
REAL :: residue = 0.
CALL READ_INTRO
CALL ludcmp(nx,5,5,1,0,0)
CALL SET_IC
pflame=0.
DO i=1,nx
pflame=pflame+y1(i,fctrl_species)*hx
ENDDO
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
t_now=0.
t_uf=0.
DO
IF (t_now.ge.tf) EXIT
residue = residual(yold, y1)
IF (converged(yold, y1)) 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(*,*)
ENDIF
ncyc=ncyc+1
t_uf=t_uf+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
! 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
ENDIF
ENDDO
! CALL write_sd
CALL write_sd
CALL write_pre
CALL save_final_field
@ -204,41 +78,40 @@
WRITE(*,*) 'Fin.'
WRITE(*,*)
END SUBROUTINE solve_
END SUBROUTINE solve
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
REAL, DIMENSION(2,nx) :: ux,dux,d2ux
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) ! Yr
ux(2,i)=y1(i,2) ! T
ux(1,i)=y1(i,1,1) ! Yr
IF (ux(1,i).gt.1.) ux(1,i)=1.
ENDDO
nd=2
CALL dfnonp(nx,hx,ux(:,:),dux(:,:),nd,1)
CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),nd,1)
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)
theta=ux(2,i)
c=1.-yr
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
! wrate=min_wr
@ -246,22 +119,22 @@
! ((exp((c-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr
! ENDIF
IF (theta.le.c_ref) THEN
IF (c.le.c_ref) THEN
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)))
ENDIF
dely=ABS(dux(1,i))
sdr=wrate/dely
sdd=-dm(i)*d2ux(1,i)/dely
sdd=-diff*d2ux(1,i)/dely
IF (dely.eq.0.) THEN
sdr=0.; sdd=0.
ENDIF
sd=sdr+sdd
uu=u(1)
uu=u(1,1,1)
onelw=(-dux(1,i))/c
dd=-dm(i)*d2ux(1,i)
dd=-diff*d2ux(1,i)
if (c.eq.0.) onelw=0.
sav(1,i)=sav(1,i)+yr
@ -282,63 +155,48 @@
END SUBROUTINE write_sd
SUBROUTINE save_final_field
INTEGER :: i
INTEGER :: i,j,k
OPEN (305,FILE='sfield.bin',form='unformatted',status='unknown')
DO i=1,nx
WRITE (305) y1(i,1)
WRITE (305) y1(i,1,1)
ENDDO
CLOSE (305)
OPEN (305,FILE='sfield.dat')
DO i=1,nx
WRITE (305,*) (/ i*hx, y1(i,:) /)
WRITE (305,'(e30.20)') y1(i,1,1)
ENDDO
CLOSE (305)
END SUBROUTINE save_final_field
SUBROUTINE write_pre
REAL :: theta,yr,c,dy,maxdy=0.,del_f
REAL :: S_L=0.,wrate,wrate1,wrate2
REAL :: ya,yx
REAL :: yr,c,dy,maxdy=0.,del_f
REAL :: S_L=0.,wrate
INTEGER :: i
REAL, DIMENSION(1,nx) :: ux, dux
if ( reaction_type == "onestep" ) then
DO i=1,nx
yr=y1(i,1)
theta=y1(i,2)
yr=y1(i,1,1)
c=1.-yr
ux(1,i)=yr
ux(1,i)=theta
wrate=rate_1step(yr, theta)
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
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
CALL dfnonp(nx,hx,ux(1,:),dux(1,:),1,1)
@ -348,338 +206,253 @@
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)
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
INTEGER :: i
DO i = 1, nsp
y1(1,i) = inletbc(i)
END DO
y1(1,1,1)=1.-ctmp
END SUBROUTINE SET_BC
SUBROUTINE BASE_FLAME_PROFILE(x)
SUBROUTINE SET_IC
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)
! initialize velocity field
u=u0;
u=u0; si=INT(nx*(1.-tar_lo))
! initial flame center (grid index)
si=INT(nx*(1.-tar_lo))
! initialize Yr field
DO i=1,nx
IF(i< nx-(si+ifl/2)) THEN
x(i)=0.+ctmp
xi=0.+ctmp
ELSE IF(i> nx-(si-ifl/2)) THEN
x(i)=1.
xi=1.
ELSE
x(i)=0.5+REAL(i-nx+si)/REAL(ifl)
xi=0.5+REAL(i-nx+si)/REAL(ifl)
ENDIF
END DO
! initialize species field
! max_ysum = maxval(y1(:,1) + y1(:,2))
! y1(:,1) = y1(:,1) / maxval(y1(:,1) + y1(:,2))
! 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))
! y1(i,1,1)=(1.-xi)*1.
y1(i,1,1)=(1.-xi) ! reactant mass fraction
ENDDO
END SUBROUTINE BASE_FLAME_PROFILE
pflame=0.
DO i=1,nx
pflame=pflame+y1(i,1,1)*hx
ENDDO
SUBROUTINE SET_IC_ONESTEP
INTEGER :: i
REAL :: xi
REAL :: x(nx)
END SUBROUTINE SET_IC
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
y1(i,2)=(xi) ! reactant mass fraction
ENDDO
END SUBROUTINE SET_IC_ONESTEP
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
SUBROUTINE SET_IC_TWOSTEP
INTEGER :: i
REAL :: xi
REAL :: x(nx)
CLOSE(itape)
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
xi = x(i)
! refwr=pre*1.*exp(-ac/(1.+bc*c_cut))
! 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) * &
exp (-(beta1*(1. - xi))/(1. - hrp*(1. - xi)))
y1(i,3) = xi
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(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.
ALLOCATE(dm(nx),STAT=ierr) ; dm=diff
if ( reaction_type == "onestep" ) then
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
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
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
INTEGER :: i,j,k
REAL :: at,bt
CALL calc_diff(ri)
CALL rhs(ri,f)
CALL fns(ri,f)
IF(istage<5) THEN
at=a(istage)*dt
bt=(b(istage)-a(istage))*dt
r1=r1+at*f
r2=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)+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
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
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)
REAL, INTENT(IN),DIMENSION(:,:) :: r1
REAL, INTENT(OUT),DIMENSION(:,:) :: f
REAL, DIMENSION(3,nx) :: ux, dux, d2ux
INTEGER :: i
REAL :: Ly,Dy,Lt,Dt
REAL :: wrate1, wrate2
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 i=1,nx
ux(1,i)=r1(i,1) ! Y
ux(2,i)=r1(i,2) ! Y
ux(3,i)=r1(i,3) ! T
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
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
wrate1=rate1_2step(ux(1,i), ux(2,i), ux(3,i))
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.
Dy=Ly
f(nx,1,1)=f(nx,1,1)-Dy
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