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

Compare commits

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

7 changed files with 55 additions and 186 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

14
code/makefile
View file

@ -1,24 +1,18 @@
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
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
touch ysolve.mod
Compact.o compact.mod : Compact.f90 compact.mod: Compact.f90
${compiler} -c ${flags} Compact.f90 ${compiler} -c ${flags} Compact.f90
touch compact.mod
test: ex test: ex
sh test.sh sh test.sh

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

@ -16,6 +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
rvis 1.0
cstar 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

135
code/ysolve.f90
View file

@ -4,14 +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 :: rvis,cstar,sigw
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,y1,y2,yf REAL, DIMENSION(:,:,:), ALLOCATABLE :: u,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., &
@ -25,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
@ -116,7 +84,6 @@
REAL :: c,yr,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(1,nx) :: ux,dux,d2ux REAL, DIMENSION(1,nx) :: ux,dux,d2ux
REAL, DIMENSION(1,nx) :: duy
REAL, DIMENSION(10,nx) :: sav REAL, DIMENSION(10,nx) :: sav
sav=0. sav=0.
@ -134,15 +101,11 @@
IF (ux(1,i).gt.1.) ux(1,i)=1. IF (ux(1,i).gt.1.) ux(1,i)=1.
ENDDO ENDDO
dm = (/ (diff * diffusivity_sutherland(1.-y1(i,1,1)), i=1,nx) /)
nd=1 nd=1
CALL dfnonp(nx,hx,ux(1,:),dux(1,:),nd,1) CALL dfnonp(nx,hx,ux(1,:),dux(1,:),nd,1)
CALL dfnonp(nx,hx,dm,duy(1,:),nd,1)
nd=1 nd=1
CALL d2fnonp(nx,hx,ux(1,:),d2ux(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)
c=1.-yr c=1.-yr
@ -164,14 +127,14 @@
dely=ABS(dux(1,i)) dely=ABS(dux(1,i))
sdr=wrate/dely sdr=wrate/dely
sdd=-(dm(i)*d2ux(1,i)+duy(1,i)*dux(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,1,1) uu=u(1,1,1)
onelw=(-dux(1,i))/c onelw=(-dux(1,i))/c
dd=-(dm(i)*d2ux(1,i)+duy(1,i)*dux(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
@ -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,rvis
WRITE(otape,*) cdum,rvis
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)
@ -364,7 +318,6 @@
ALLOCATE(yf(nx,1,1),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
@ -411,7 +364,6 @@
REAL, INTENT(IN),DIMENSION(:,:,:) :: r1 REAL, INTENT(IN),DIMENSION(:,:,:) :: r1
REAL, INTENT(OUT),DIMENSION(:,:,:) :: f REAL, INTENT(OUT),DIMENSION(:,:,:) :: f
REAL, DIMENSION(3,nx) :: ux,dux,d2ux REAL, DIMENSION(3,nx) :: ux,dux,d2ux
REAL, DIMENSION(1,nx) :: duy
INTEGER :: i,j,k INTEGER :: i,j,k
REAL :: wrate,Ly,Dy REAL :: wrate,Ly,Dy
@ -424,10 +376,7 @@
ux(3,i)=r1(i,j,k) ! Y ux(3,i)=r1(i,j,k) ! Y
ENDDO ENDDO
dm = (/ (diff * diffusivity_sutherland(1.-r1(i,1,1)), i=1,nx) /)
CALL dfnonp(nx,hx,ux(1:3,:),dux(1:3,:),3,1) CALL dfnonp(nx,hx,ux(1:3,:),dux(1:3,:),3,1)
CALL dfnonp(nx,hx,dm,duy(1,:),1,1)
CALL d2fnonp(nx,hx,ux(3:3,:),d2ux(1,:),1,1) CALL d2fnonp(nx,hx,ux(3:3,:),d2ux(1,:),1,1)
DO i=1,nx DO i=1,nx
@ -446,14 +395,57 @@
! -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,j,k)=-0.5*( dux(1,i) + ux(2,i)*dux(3,i) + & f(i,j,k)=-0.5*( dux(1,i) + ux(2,i)*dux(3,i) + &
ux(3,i)*dux(2,i) ) & ux(3,i)*dux(2,i) ) &
+ dm(i)*d2ux(1,i) + duy(1,i)*dux(3,i) - wrate + diff*d2ux(1,i) - wrate
IF (i.eq.nx) THEN IF (i.eq.nx) THEN
Ly=ux(2,nx)*dux(3,nx) ! Ly = u*dYr/dx Ly=ux(2,nx)*dux(3,nx) ! Ly = u*dYr/dx
f(nx,1,1)=-wrate f(nx,1,1)=-wrate
ENDIF ENDIF
ENDDO 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 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,1)=0. f(1,1,1)=0.
@ -462,34 +454,5 @@
f(nx,1,1)=f(nx,1,1)-Dy f(nx,1,1)=f(nx,1,1)-Dy
END SUBROUTINE fns END SUBROUTINE fns
!------------------------------------------------------------------------
REAL FUNCTION sigmoid(x)
REAL, INTENT(IN) :: x
sigmoid = 1. / (1 + exp(-x))
END FUNCTION sigmoid
!------------------------------------------------------------------------
REAL FUNCTION diffusivity_sigmoid(c)
REAL, INTENT(IN) :: c
diffusivity_sigmoid = 1. + (rvis - 1.) * sigmoid(sigw*(c-cstar))
END FUNCTION diffusivity_sigmoid
!------------------------------------------------------------------------
REAL FUNCTION diffusivity_T076(c)
REAL, INTENT(IN) :: c
diffusivity_T076 = ((1.0 + bc * c) ** 0.76)
END FUNCTION diffusivity_T076
!------------------------------------------------------------------------
REAL FUNCTION diffusivity_sutherland(c)
REAL, INTENT(IN) :: c
REAL :: theta, As, Ts, T0, T1
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 END MODULE ysolve