diff --git a/code/m_chemistry.f90 b/code/m_chemistry.f90 index 26845cd..67024b2 100644 --- a/code/m_chemistry.f90 +++ b/code/m_chemistry.f90 @@ -14,6 +14,18 @@ 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 @@ -22,6 +34,7 @@ subroutine init_chemistry beta1_twostep = beta1 hrp_twostep = hrp else + WRITE(*,*) 'ERROR, UNDEFINED REACTION TYPE ', reaction_type stop end if diff --git a/code/ysolve.f90 b/code/ysolve.f90 index 8ab3149..4bdb072 100644 --- a/code/ysolve.f90 +++ b/code/ysolve.f90 @@ -6,6 +6,7 @@ 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 @@ -121,7 +122,7 @@ pflame=0. DO i=1,nx pflame=pflame+y1(i,1)*hx - uxt(i)=y1(i,3) + uxt(i)= 1. - y1(i,1) ENDDO CALL dfnonp(nx,hx,uxt,duxt,1,1) @@ -287,14 +288,17 @@ END SUBROUTINE write_pre SUBROUTINE SET_BC - y1(1,1)=1. - y1(1,2)=0. - y1(1,3)=0. + 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) @@ -309,22 +313,37 @@ ! initialize Yr field DO i=1,nx IF(i< nx-(si+ifl/2)) THEN - xi=0.+ctmp + x(i)=0.+ctmp ELSE IF(i> nx-(si-ifl/2)) THEN - xi=1. + x(i)=1. ELSE - xi=0.5+REAL(i-nx+si)/REAL(ifl) + x(i)=0.5+REAL(i-nx+si)/REAL(ifl) ENDIF -! y1(i,1,1)=(1.-xi)*1. - y1(i,1)=(1.-xi) ! reactant mass fraction + END DO - y1(i,2) = (1./2.) * (lambda1/lambda2) * y1(i,1) * & - exp (-(beta1*(1. - xi))/(1. - hrp*(1. - xi))) + if ( reaction_type == "onestep" ) then + DO i=1,nx + xi = x(i) - y1(i,3) = xi + 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 - ENDDO ! initialize species field ! max_ysum = maxval(y1(:,1) + y1(:,2)) @@ -350,14 +369,29 @@ CALL init_chemistry ALLOCATE( u(nx),STAT=ierr) ; u=0. - ALLOCATE( y1(nx,3),STAT=ierr) ; y1=0. - ALLOCATE( y2(nx,3),STAT=ierr) ; y2=0. - ALLOCATE( yf(nx,3),STAT=ierr) ; yf=0. - ALLOCATE(yold(nx,3),STAT=ierr) ; yold=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 @@ -365,6 +399,7 @@ SUBROUTINE finalize_solver DEALLOCATE( u) + DEALLOCATE(inletbc) DEALLOCATE(y1) DEALLOCATE(y2) DEALLOCATE(yf)