516 lines
13 KiB
Fortran
516 lines
13 KiB
Fortran
MODULE ysolve
|
|
USE m_compact
|
|
USE m_parameters
|
|
USE m_chemistry
|
|
|
|
IMPLICIT NONE
|
|
|
|
REAL, DIMENSION(:), ALLOCATABLE :: u
|
|
REAL, DIMENSION(:), ALLOCATABLE :: inletbc
|
|
REAL, DIMENSION(:,:), ALLOCATABLE :: y1,y2,yf,yold
|
|
REAL, DIMENSION(:), ALLOCATABLE :: uxt,duxt
|
|
REAL, DIMENSION(:), ALLOCATABLE :: dm
|
|
|
|
CONTAINS
|
|
!------------------------------------------------------------------------
|
|
SUBROUTINE parse
|
|
CHARACTER(100) :: num1char
|
|
|
|
!First, make sure the right number of inputs have been provided
|
|
IF(COMMAND_ARGUMENT_COUNT().EQ.1)THEN
|
|
CALL GET_COMMAND_ARGUMENT(1,num1char) !first, read in the two values
|
|
ELSE IF(COMMAND_ARGUMENT_COUNT().EQ.0)THEN
|
|
WRITE(*,*)'ERROR, NO COMMAND-LINE ARGUMENT'
|
|
num1char = "-h"
|
|
ELSE
|
|
WRITE(*,*)'ERROR, TOO MANY COMMAND-LINE ARGUMENTS(MORE THAN ONE). STOPPING'
|
|
STOP
|
|
ENDIF
|
|
|
|
|
|
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
|
|
|
|
ELSE
|
|
|
|
END IF
|
|
|
|
|
|
|
|
|
|
IF(read_stdin)THEN
|
|
WRITE(*,*) "Read from STDIN"
|
|
ELSE
|
|
itape_name = num1char
|
|
WRITE(*,*) "Read from " // 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
|
|
|
|
INTEGER :: i,j,k,savenum
|
|
REAL :: pflame,pflold,delf=0.
|
|
REAL :: residue = 0.
|
|
|
|
CALL SET_IC
|
|
|
|
t_now=0.
|
|
t_uf=0.
|
|
|
|
DO
|
|
|
|
IF (t_now.ge.tf) EXIT
|
|
|
|
residue = residual(yold, y1)
|
|
|
|
IF (converged(yold, y1)) EXIT
|
|
|
|
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(fns)
|
|
|
|
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
|
|
ENDDO
|
|
|
|
CALL write_sd
|
|
CALL write_pre
|
|
CALL save_final_field
|
|
|
|
WRITE(*,*)
|
|
WRITE(*,*) 'Fin.'
|
|
WRITE(*,*)
|
|
|
|
END SUBROUTINE solve
|
|
|
|
SUBROUTINE write_sd
|
|
REAL :: c,yr,theta,wrate,dely,sdr,sdd,sd,uu,onelw,dd
|
|
INTEGER :: i,j,k,nd
|
|
REAL, DIMENSION(2,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
|
|
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)
|
|
|
|
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
|
|
|
|
! 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 (theta.le.c_ref) THEN
|
|
wrate=min_wr
|
|
IF (theta.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(theta-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
|
|
IF (dely.eq.0.) THEN
|
|
sdr=0.; sdd=0.
|
|
ENDIF
|
|
sd=sdr+sdd
|
|
uu=u(1)
|
|
onelw=(-dux(1,i))/c
|
|
dd=-dm(i)*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)
|
|
ENDDO
|
|
CLOSE (305)
|
|
|
|
OPEN (305,FILE='sfield.dat')
|
|
DO i=1,nx
|
|
WRITE (305,'(4e30.20)') i*hx, y1(i,1), y1(i,2), y1(i,3)
|
|
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
|
|
INTEGER :: i
|
|
REAL, DIMENSION(1,nx) :: ux, dux
|
|
|
|
DO i=1,nx
|
|
yr=y1(i,1)
|
|
|
|
ya=y1(i,1)
|
|
yx=y1(i,2)
|
|
theta=y1(i,3)
|
|
|
|
ux(1,i)=theta
|
|
|
|
wrate=rate_1step(yr, theta)
|
|
|
|
wrate1=rate1_2step(ya, yx, theta)
|
|
|
|
wrate2=rate2_2step(yx, theta)
|
|
|
|
S_L=S_L+wrate1*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)
|
|
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
|
|
END SUBROUTINE SET_BC
|
|
|
|
SUBROUTINE SET_IC
|
|
INTEGER :: i, ifl, si
|
|
REAL :: xi
|
|
REAL :: x(nx)
|
|
REAL :: max_ysum
|
|
|
|
! initial flame thickness (0.2 pi n grids)
|
|
ifl=INT(2.0*pi*0.1/hx)
|
|
|
|
! initialize velocity field
|
|
u=u0;
|
|
|
|
! 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
|
|
ELSE IF(i> nx-(si-ifl/2)) THEN
|
|
x(i)=1.
|
|
ELSE
|
|
x(i)=0.5+REAL(i-nx+si)/REAL(ifl)
|
|
ENDIF
|
|
|
|
END DO
|
|
|
|
if ( reaction_type == "onestep" ) then
|
|
DO i=1,nx
|
|
xi = x(i)
|
|
|
|
y1(i,1)=(1.-xi) ! reactant mass fraction
|
|
ENDDO
|
|
else if ( reaction_type == "twostep" ) then
|
|
DO i=1,nx
|
|
xi = x(i)
|
|
|
|
y1(i,1)=(1.-xi) ! reactant mass fraction
|
|
|
|
y1(i,2) = (1./2.) * (lambda1/lambda2) * y1(i,1) * &
|
|
exp (-(beta1*(1. - xi))/(1. - hrp*(1. - xi)))
|
|
|
|
y1(i,3) = xi
|
|
ENDDO
|
|
else
|
|
WRITE(*,*) 'ERROR, UNDEFINED REACTION TYPE ', reaction_type
|
|
stop
|
|
end if
|
|
|
|
|
|
! 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))
|
|
|
|
! calculate flame position
|
|
pflame=0.
|
|
DO i=1,nx
|
|
pflame=pflame+y1(i,fctrl_species)*hx
|
|
ENDDO
|
|
|
|
END SUBROUTINE SET_IC
|
|
|
|
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(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
|
|
|
|
|
|
|
|
|
|
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)
|
|
|
|
INTERFACE
|
|
SUBROUTINE rhs(r1,f)
|
|
REAL, INTENT(IN), DIMENSION(:,:) :: r1
|
|
REAL, INTENT(OUT), DIMENSION(:,:) :: f
|
|
END SUBROUTINE rhs
|
|
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 rhs(ri,f)
|
|
|
|
IF(istage<5) THEN
|
|
at=a(istage)*dt
|
|
bt=(b(istage)-a(istage))*dt
|
|
r1=r1+at*f
|
|
r2=r1+bt*f
|
|
ELSE
|
|
bt=b(istage)*dt
|
|
r1=r1+bt*f
|
|
ENDIF
|
|
END SUBROUTINE substep
|
|
|
|
END SUBROUTINE RK4
|
|
!------------------------------------------------------------------------
|
|
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,Lt,Dt
|
|
REAL :: wrate1, wrate2
|
|
|
|
! 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
|
|
ENDDO
|
|
|
|
CALL dfnonp(nx,hx,ux(:,:),dux(:,:),3,1)
|
|
|
|
CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),3,1)
|
|
|
|
DO i=1,nx
|
|
wrate=rate_1step(ux(1,i), ux(3,i))
|
|
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.
|
|
|
|
END SUBROUTINE fns
|
|
|
|
!------------------------------------------------------------------------
|
|
|
|
END MODULE ysolve
|