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

Compare commits

base: research:nonunityle
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. "nonunityle" and "master" have entirely different histories.
nonunityle ... master

7 changed files with 130 additions and 220 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

16
code/makefile
View file

@ -1,25 +1,17 @@
flags = -Wall -O3 -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8 -march=native flags = -Wall -O3 -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 ysolve.o Compact.o ex : test.o ysolve.mod compact.mod
${compiler} -o ex test.o ysolve.o Compact.o ${compiler} -o ex test.o ysolve.o Compact.o
test.o : test.f90 ysolve.mod compact.mod test.o : test.f90 ysolve.mod
${compiler} -c ${flags} test.f90 ${compiler} -c ${flags} test.f90
ysolve.o ysolve.mod : ysolve.f90 compact.mod ysolve.mod: ysolve.f90 compact.mod
${compiler} -c ${flags} ysolve.f90 ${compiler} -c ${flags} ysolve.f90
Compact.o compact.mod : Compact.f90 compact.mod: Compact.f90
${compiler} -c ${flags} Compact.f90 ${compiler} -c ${flags} Compact.f90
test: ex test: ex

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

@ -16,5 +16,4 @@
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

83
code/sample/itape.py
View file

@ -1,83 +0,0 @@
#!/bin/python
import os, sys
import argparse
parser = argparse.ArgumentParser()
file_parser = argparse.ArgumentParser()
param_dict = {
"nx" :512,
"l_0" :4.,
"int_pr" :400,
"tar_lo" :0.60,
"dt" :0.001,
"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,
"rvis" :1.0,
"cstar" :1.0,
"sigw" :100.0,
}
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)
itape = ''' nx {0[nx]}
l_0 {0[l_0]}
int_pr {0[int_pr]}
tar_lo {0[tar_lo]}
dt {0[dt]}
sc {0[sc]}
vis {0[vis]}
pre {0[pre]}
ac {0[ac]}
bc {0[bc]}
u0 {0[u0]}
tf {0[tf]}
dt_uf {0[dt_uf]}
ctmp {0[ctmp]}
c_cut {0[c_cut]}
c_ref {0[c_ref]}
min_wr {0[min_wr]}
prof_wr {0[prof_wr]}
rvis {0[rvis]}
cstar {0[cstar]}
sigw {0[sigw]} '''
print itape.format( params )
# 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")

1
code/test.f90
View file

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

2
code/test.sh
View file

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

243
code/ysolve.f90
View file

@ -4,15 +4,12 @@
PRIVATE PRIVATE
REAL, PARAMETER :: pi=3.14159265358979323846 REAL, PARAMETER :: pi=3.14159265358979323846
REAL :: hx,dt,vis,sc,diff,pre,ac,bc,tf,t_now,t_uf,dt_uf REAL :: hx,dt,vis,sc,diff,pre,ac,bc,tf,t_now,t_uf,dt_uf
REAL :: lewis
REAL :: c_cut,c_ref,refwr,minf,tar_lo,u0,ctmp,l_0,lo_flm=0. 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 :: er_lof=0., erdot=0.,min_wr,prof_wr
REAL :: pflame,pflold,oldu REAL :: pflame,pflold,oldu
INTEGER :: ncyc=0,int_pr,nx INTEGER :: ncyc=0,int_pr,nx
REAL, DIMENSION(:), ALLOCATABLE :: u REAL, DIMENSION(:,:,:), ALLOCATABLE :: u,y1,y2,yf
REAL, DIMENSION(:,:), ALLOCATABLE :: y1,y2,yf
REAL, DIMENSION(:), ALLOCATABLE :: uxt,duxt REAL, DIMENSION(:), ALLOCATABLE :: uxt,duxt
REAL, DIMENSION(:), ALLOCATABLE :: dm
INTEGER :: istage INTEGER :: istage
REAL, DIMENSION(5) :: a=(/ 970286171893./4311952581923., & REAL, DIMENSION(5) :: a=(/ 970286171893./4311952581923., &
@ -26,38 +23,8 @@
2114624349019./3568978502595., & 2114624349019./3568978502595., &
5198255086312./14908931495163. /) 5198255086312./14908931495163. /)
LOGICAL :: read_itape, read_stdin PUBLIC :: solve
CHARACTER(100) :: itape_name
PUBLIC :: solve, parse
CONTAINS CONTAINS
!------------------------------------------------------------------------
SUBROUTINE parse
CHARACTER(100) :: num1char
!First, make sure the right number of inputs have been provided
IF(COMMAND_ARGUMENT_COUNT().EQ.0)THEN
! read from itape
read_itape = .true.
ELSE IF(COMMAND_ARGUMENT_COUNT().EQ.1)THEN
CALL GET_COMMAND_ARGUMENT(1,num1char) !first, read in the two values
ELSE
WRITE(*,*)'ERROR, TOO MANY COMMAND-LINE ARGUMENTS(MORE THAN ONE). STOPPING'
STOP
ENDIF
read_stdin = (num1char=="-")
IF(read_itape)THEN
itape_name = "itape"
WRITE(*,*) "Read from itape"
ELSE IF(read_stdin)THEN
WRITE(*,*) "Read from STDIN"
ELSE
itape_name = num1char
WRITE(*,*) "Read from " // num1char
END IF
END SUBROUTINE parse
!------------------------------------------------------------------------ !------------------------------------------------------------------------
SUBROUTINE solve SUBROUTINE solve
INTEGER :: i,j,k,savenum INTEGER :: i,j,k,savenum
@ -79,17 +46,17 @@
pflold=pflame pflold=pflame
pflame=0. pflame=0.
DO i=1,nx DO i=1,nx
pflame=pflame+y1(i,1)*hx pflame=pflame+y1(i,1,1)*hx
uxt(i)=y1(i,1) uxt(i)=y1(i,1,1)
ENDDO ENDDO
CALL dfnonp(nx,hx,uxt,duxt,1,1) CALL dfnonp(nx,hx,uxt,duxt,1,1)
delf=1./MAXVAL(ABS(duxt)) delf=1./MAXVAL(ABS(duxt))
oldu=u(1) oldu=u(1,1,1)
u=u+0.5*(hx*REAL(nx)*tar_lo-pflame)+0.5*(pflold-pflame) u=u+0.5*(hx*REAL(nx)*tar_lo-pflame)+0.5*(pflold-pflame)
t_uf=0. t_uf=0.
WRITE(*,'(a3,f8.3,a10,f6.3,a10,f6.3,a10,f7.4,a10,f7.4,a8,f7.4)') & 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, & ' 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 ' // Old_U:',oldu,' // New_U:',u(1,1,1),' // L_f:',delf
! 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(*,'(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(*,*) ! WRITE(*,*)
@ -114,9 +81,9 @@
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.
@ -130,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
@ -153,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
@ -193,39 +159,38 @@
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,'(2e30.20)') y1(i,1), y1(i,2) 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 REAL :: S_L=0.,wrate
INTEGER :: i INTEGER :: i
REAL, DIMENSION(1,nx) :: ux, dux REAL, DIMENSION(1,nx) :: ux, dux
DO i=1,nx DO i=1,nx
yr=y1(i,1) yr=y1(i,1,1)
theta=y1(i,2)
c=1.-yr c=1.-yr
ux(1,i)=yr ux(1,i)=yr
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
! IF (c.gt.c_ref) wrate= & ! IF (c.gt.c_ref) wrate= &
! ((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
@ -242,7 +207,7 @@
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(*,*)
@ -250,7 +215,7 @@
SUBROUTINE SET_BC SUBROUTINE SET_BC
y1(1,1)=1.-ctmp y1(1,1,1)=1.-ctmp
END SUBROUTINE SET_BC END SUBROUTINE SET_BC
SUBROUTINE SET_IC SUBROUTINE SET_IC
@ -271,14 +236,12 @@
ENDIF ENDIF
! y1(i,1,1)=(1.-xi)*1. ! y1(i,1,1)=(1.-xi)*1.
y1(i,1)=(1.-xi) ! reactant mass fraction y1(i,1,1)=(1.-xi) ! reactant mass fraction
ENDDO ENDDO
y1(:,2) = 1.0 - y1(:,1)
pflame=0. pflame=0.
DO i=1,nx DO i=1,nx
pflame=pflame+y1(i,1)*hx pflame=pflame+y1(i,1,1)*hx
ENDDO ENDDO
END SUBROUTINE SET_IC END SUBROUTINE SET_IC
@ -286,14 +249,9 @@
SUBROUTINE READ_INTRO SUBROUTINE READ_INTRO
CHARACTER(LEN=8) :: cdum CHARACTER(LEN=8) :: cdum
INTEGER :: itape=300,otape=301,ierr INTEGER :: itape=300,otape=301,ierr
INTEGER :: i
IF (read_stdin) THEN
itape=5
ELSE
OPEN(itape,FILE=itape_name)
END IF
OPEN(itape,FILE='itape')
OPEN(otape,FILE='otape') OPEN(otape,FILE='otape')
READ(itape,*) cdum,nx READ(itape,*) cdum,nx
WRITE(otape,*) cdum,nx WRITE(otape,*) cdum,nx
@ -333,12 +291,8 @@
WRITE(otape,*) cdum,min_wr WRITE(otape,*) cdum,min_wr
READ(itape,*) cdum,prof_wr READ(itape,*) cdum,prof_wr
WRITE(otape,*) cdum,prof_wr WRITE(otape,*) cdum,prof_wr
READ(itape,*) cdum,lewis
WRITE(otape,*) cdum,lewis
IF (.not.read_stdin) THEN CLOSE(itape)
CLOSE(itape)
END IF
! hx=l_0*pi/REAL(nx) ! hx=l_0*pi/REAL(nx)
hx=l_0*pi/REAL(nx-1) hx=l_0*pi/REAL(nx-1)
@ -358,13 +312,12 @@
refwr=pre*1.*exp(-ac/(1.+bc*c_ref)) refwr=pre*1.*exp(-ac/(1.+bc*c_ref))
l_0=l_0*pi l_0=l_0*pi
ALLOCATE( u(nx),STAT=ierr) ; u=0. ALLOCATE(u(nx,1,1),STAT=ierr) ; u=0.
ALLOCATE(y1(nx,2),STAT=ierr) ; y1=0. ALLOCATE(y1(nx,1,1),STAT=ierr) ; y1=0.
ALLOCATE(y2(nx,2),STAT=ierr) ; y2=0. ALLOCATE(y2(nx,1,1),STAT=ierr) ; y2=0.
ALLOCATE(yf(nx,2),STAT=ierr) ; yf=0. ALLOCATE(yf(nx,1,1),STAT=ierr) ; yf=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
END SUBROUTINE READ_INTRO END SUBROUTINE READ_INTRO
@ -377,8 +330,8 @@
END SUBROUTINE RK4 END SUBROUTINE RK4
!------------------------------------------------------------------------ !------------------------------------------------------------------------
SUBROUTINE substep(ri,r1,r2,f) SUBROUTINE substep(ri,r1,r2,f)
REAL, INTENT(INOUT),DIMENSION(:,:) :: ri,r1,r2 REAL, INTENT(INOUT),DIMENSION(:,:,:) :: ri,r1,r2
REAL, INTENT(OUT),DIMENSION(:,:) :: f REAL, INTENT(OUT),DIMENSION(:,:,:) :: f
INTEGER :: i,j,k INTEGER :: i,j,k
REAL :: at,bt REAL :: at,bt
@ -387,66 +340,118 @@
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
!------------------------------------------------------------------------ !------------------------------------------------------------------------
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(2,nx) :: ux, dux, d2ux REAL, DIMENSION(3,nx) :: ux,dux,d2ux
INTEGER :: i,j,k INTEGER :: i,j,k
REAL :: wrate,Ly,Dy,Lt,Dt REAL :: wrate,Ly,Dy
DO k=1,1 !nz
! x-direction ! x-direction
DO j=1,1 !ny
DO i=1,nx DO i=1,nx
ux(1,i)=r1(i,1) ! Y ux(1,i)=u(i,j,k)*r1(i,j,k) ! u*Y
ux(2,i)=r1(i,2) ! T ux(2,i)=u(i,j,k) ! u
ux(3,i)=r1(i,j,k) ! Y
ENDDO ENDDO
CALL dfnonp(nx,hx,ux(:,:),dux(:,:),2,1) CALL dfnonp(nx,hx,ux(1:3,:),dux(1:3,:),3,1)
CALL d2fnonp(nx,hx,ux(3:3,:),d2ux(1,:),1,1)
CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),2,1)
DO i=1,nx DO i=1,nx
wrate=pre*ux(1,i)*exp(-ac/(1.+bc*(ux(2,i)))) !wrate wrate=pre*ux(3,i)*exp(-ac/(1.+bc*(1.-ux(3,i)))) !wrate
! IF ((1.-ux(3,i)).le.c_cut) THEN
IF ((ux(2,i)).le.c_ref) 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 wrate=min_wr
IF ((ux(2,i)).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(ux(2,i)-c_ref))+ & 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))) min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
ENDIF ENDIF
! -0.5*( d(u*Y)/dx + u*dY/dx + Y*du/dx ) + D*d2Y/d2x ! -0.5*( d(u*Y)/dx + u*dY/dx + Y*du/dx ) + D*d2Y/d2x
f(i,1) = - ( u(i)*dux(1,i) ) + diff * d2ux(1,i) - wrate f(i,j,k)=-0.5*( dux(1,i) + ux(2,i)*dux(3,i) + &
ux(3,i)*dux(2,i) ) &
! - u*dT/dx + D*d2T/d2x + diff*d2ux(1,i) - wrate
f(i,2) = - ( u(i)*dux(2,i) ) + (lewis * diff) * d2ux(2,i) + wrate
IF (i.eq.nx) THEN IF (i.eq.nx) THEN
Ly=u(nx)*dux(1,nx) ! Ly = u*dYr/dx Ly=ux(2,nx)*dux(3,nx) ! Ly = u*dYr/dx
f(nx,1)=-wrate f(nx,1,1)=-wrate
Lt=u(nx)*dux(2,nx) ! Ly = u*dYr/dx
f(nx,2)=wrate
ENDIF ENDIF
ENDDO 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
! 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
! Boundary conditionS ! Boundary conditionS
f(1,1)=0. f(1,1,1)=0.
Dy=Ly Dy=Ly
f(nx,1)=f(nx,1)-Dy f(nx,1,1)=f(nx,1,1)-Dy
f(1,2)=0.
Dt=Lt
f(nx,2)=f(nx,2)-Dt
END SUBROUTINE fns END SUBROUTINE fns
!------------------------------------------------------------------------ !------------------------------------------------------------------------