diff --git a/code/ysolve.f90 b/code/ysolve.f90 index d967215..37d8ef9 100644 --- a/code/ysolve.f90 +++ b/code/ysolve.f90 @@ -90,22 +90,25 @@ SUBROUTINE solve IF (reaction_type == "onestep") THEN - CALL solve_ (fonestep) + CALL solve_ (fonestep, SET_IC_ONESTEP) ELSE IF (reaction_type == "twostep") THEN - CALL solve_ (fns) + CALL solve_ (fns, SET_IC_TWOSTEP) ELSE WRITE(*,*)'ERROR, UNSUPPORTED REACTION ', trim(reaction_type), '. STOPPING' STOP END IF END SUBROUTINE solve - SUBROUTINE solve_ (rhs) + SUBROUTINE solve_ (rhs, SET_IC) 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 END INTERFACE INTEGER :: i,j,k,savenum @@ -114,6 +117,12 @@ CALL SET_IC + pflame=0. + DO i=1,nx + pflame=pflame+y1(i,fctrl_species)*hx + ENDDO + + t_now=0. t_uf=0. @@ -331,7 +340,7 @@ END DO END SUBROUTINE SET_BC - SUBROUTINE SET_IC + SUBROUTINE BASE_FLAME_PROFILE(x) INTEGER :: i, ifl, si REAL :: xi REAL :: x(nx) @@ -355,18 +364,39 @@ ELSE x(i)=0.5+REAL(i-nx+si)/REAL(ifl) ENDIF - END DO - if ( reaction_type == "onestep" ) then + ! 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)) + + END SUBROUTINE BASE_FLAME_PROFILE + + SUBROUTINE SET_IC_ONESTEP + INTEGER :: i + REAL :: xi + REAL :: x(nx) + + CALL BASE_FLAME_PROFILE(x) + DO i=1,nx xi = x(i) y1(i,1)=(1.-xi) ! reactant mass fraction y1(i,2)=(xi) ! reactant mass fraction - ENDDO - else if ( reaction_type == "twostep" ) then + END SUBROUTINE SET_IC_ONESTEP + + SUBROUTINE SET_IC_TWOSTEP + INTEGER :: i + REAL :: xi + REAL :: x(nx) + + CALL BASE_FLAME_PROFILE(x) + DO i=1,nx xi = x(i) @@ -377,26 +407,7 @@ 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 + END SUBROUTINE SET_IC_TWOSTEP SUBROUTINE init_solver