c =====================================================
subroutine combl()
c =====================================================
c
c Create and initialize application specific common-blocks.
c
implicit none
include 'cles.i'
include 'heairacesf6.i'
include 'cuser.i'
integer lin,lout, i,nv,m,k
double precision tmp
double precision vec(7)
double precision gm1, gp1, g2, M2,MachNo
double precision c_air, c_Sair, c_Usf6
double precision vx(7)
double precision gamma1_air,gamma1_sf6,gamma1_hemx
double precision gamma_hemx,rgas_hemx,c_hemx
double precision pwr
double precision x,th
parameter( lin = 5, lout = 6 )
nv = 7
call stdcdat()
c
c GET INIT. CONDITIONS
c read in the density of air and the ambient pressure
c (the unshocked state of air).
c Uniform abient pressure and tempurature defines the
c density of sf_6
c read in the location of the shock and the contact
c read in the mach number of the shock
open(unit=lin,status='old',form='formatted',file='init.dat')
c the shock number and location of diaphragm
read (lin, *) MachNo,dloc
c ambient temp and pressure
read (lin, *) t_air,p_air
c driver gas temp
read (lin, *) t_hemx,fraction_air
c the locations and thickness of the sides of the sf6 curtain
read (lin, *) crtnL,crtnR
read (lin, *) zicki
c shock capturing info
c read (lin, *) dwidth, minc, minjump, gradS
read (lin, *) CurvP,CurvRho,CurvS
c perturbation fequencies on the left and right sides
read (lin, *) omegaL,omegaR
c perturbation amplitudes
read (lin, *) ampL,ampR
c trun=1, for smaller, with shock intialized at x=sloc
c trun=0, is the full tube.
read (lin, *) trunc,sloc
close (lin)
c # test that the fraction of air in the driver is less then 1.
if( (fraction_air.lt.0.d0).or.(fraction_air.gt.1) ) then
print *, 'the fraction of air', fraction_air, ' not in [0,1]'
stop
end if
C----------------------------------------------------------
C CALCULATE THE STATE OF unshocked AIR
C----------------------------------------------------------
c
c the density of the unshocked air
rho_air = p_air/(rgas_air*t_air)
c the unshocked gas is at rest
u_air = 0.d0
gamma1_air = gamma_air-1.d0
c_air = dsqrt(gamma_air*rgas_air*t_air)
rhou_air = rho_air * u_air
e_air = p_air/gamma1_air + 0.5d0*rhou_air**2/rho_air
C----------------------------------------------------------
C CALCULATE THE STATE OF (unshocked) SF6
C----------------------------------------------------------
c assume the sf6 is at the same tempurature
c and pressure as the unshocked air
c compute the density of sf6
p_sf6 = p_air
t_sf6 = t_air
u_sf6 = u_air
rho_sf6=p_sf6/t_air/rgas_sf6
c # data in gas curtain:
gamma1_sf6 = gamma_sf6-1.d0
rhou_sf6 = rho_sf6 * u_sf6
e_sf6 = p_sf6/gamma1_sf6 + 0.5d0*rhou_sf6**2/rho_sf6
c
C----------------------------------------------------------
C CALCULATE THE STATE OF (Pressurized) Driver He
C -- this is actually a mixture of He and a small fraction
c -- of AIR
C----------------------------------------------------------
c # get rgas and gamma for the driver mixture
vx(1) = 1.d0
vx(6) = fraction_air
vx(7) = 0.d0
call cles_gamma(vx,nv,gamma_hemx,rgas_hemx)
c # compute the pressure needed in He to make the required shock.
c_hemx = sqrt(gamma_hemx*rgas_hemx*t_hemx)
c # driver initially at rest
u_hemx = 0.d0
pwr = 2.d0*gamma_hemx/(gamma_hemx-1.d0)
p_hemx = (gamma_hemx-1.d0)/(gamma_air+1.d0) * (c_air/c_hemx)
p_hemx = p_hemx*(MachNo-1.d0/MachNo)
p_hemx = (1.d0-p_hemx)**(-pwr)
p_hemx = p_hemx*(2.d0*gamma_air * MachNo*MachNo-gamma_air+1.d0)
p_hemx = p_hemx/(gamma_air+1.d0)
p_hemx = p_air*p_hemx
c # compute the densities of the He and AIR
rho_hemx = p_hemx/(rgas_hemx*t_hemx)
c # data in left state: - Hot(?) He
gamma1_hemx = gamma_hemx-1.d0
rhou_hemx = rho_hemx * u_hemx
e_hemx = p_hemx/gamma1_hemx + 0.5d0*rhou_hemx**2/rho_hemx
C----------------------------------------------------------
c CALCULATE THE STATE OF (shocked) AIR
C - we will only need this in the truncated case
C----------------------------------------------------------
c compute the R-H jump conditions for this shock
c _air: unshocked
c _Sair: shocked
M2 = MachNo*MachNo
gm1 = gamma_air-1.d0
gp1 = gamma_air+1.0d0
g2 = 2.d0*gamma_air/(gamma_air+1.d0)
c speed of sound in unshocked air
c_air = dsqrt(gamma_air*rgas_air*t_air)
c the shocked state of air
p_Sair = p_air*(1.0d0+g2*(M2-1.d0))
rho_Sair = rho_air*gp1*M2/(gm1*M2+2.0d0)
u_Sair = u_air + c_air*2.0d0*(M2-1.0d0)/MachNo/gp1
rhou_Sair = rho_Sair*u_Sair
e_Sair = p_Sair/gamma1_air + 0.5d0*rhou_Sair**2/rho_Sair
t_Sair = p_Sair/rho_Sair/rgas_air
c_Sair = dsqrt(gamma_air*p_Sair/rho_Sair)
c FOR THE TRUNCATED DOMAIN - set up Acoustic BCs
IF(trunc.EQ.1) THEN
c
c read in the acoustic BCs
open(unit=lin,status='old',form='formatted',file='ABCinit.dat')
read (lin, *) n_bc
do i=1, n_bc
read (lin, *) time_bc(i)
read (lin, *) rho_bc(i)
read (lin, *) u_bc(i)
read (lin, *) p_bc(i)
enddo
close(lin)
c test that the shock is in the air
if(sloc.gt.crtnL) then
print *, 'error, this shock starts in the sf6'
stop
end if
C----------------------------------------------------------
c Characteristic Boundary Conditions
C----------------------------------------------------------
call UseAcousticBC(CLES_CBC_MODE1D)
c Outflow
call SetAcousticBoundary(CLES_CBC_XDIR, CLES_CBC_LEFT,
$ CLES_CBC_INFLOW)
do i=1, 7
vec(i) = 1.5d0
enddo
call SetAcousticForceMatrix(CLES_CBC_XDIR, CLES_CBC_LEFT,
$ vec)
call SetAcousticMach(MachNo)
else if ( trunc .eq. 2 ) then
c # Load initial conditions from file
c # initialize up to sloc with the data in the file
open(unit=lin, file='initialc.dat', status='old',
$ form='formatted')
nic = 0
do while ( .TRUE. )
read(lin,*,end=1000) (qin(k,nic+1),k=1,9)
nic = nic + 1
if ( nic .eq. max_nic ) then
print *, 'Temporary storage array for initial field'
$ //' is too small'
stop
endif
enddo
1000 close(lin)
end if
c Flagging discontinuity criteria stuff
c Building a centered tanh
x = 0.0d0
th = 0.0d0
sumth = 0.0d0
sumth2 = 0.0d0
do m=-span,span, 1
x = m/dwidth
th = tanh(x)
thvec(m) = th
sumth = sumth + th
sumth2 = sumth2 + th*th
enddo
c Building a non-centered tanh
x = 0.0d0
th = 0.0d0
sumthbis = 0.0d0
sumth2bis = 0.0d0
do m=-span,span-1, 1
x = (m + 1.0d0/2.0d0)/dwidth
th = tanh(x)
thvecbis(m) = th
sumthbis = sumthbis + th
sumth2bis = sumth2bis + th*th
enddo
return
end
c =====================================================
subroutine comblgm(dimensions,shape,geometry)
c =============================+=======================
c
c Create and initialize application specific common-blocks.
c
implicit none
c
include 'cuser.i'
c Arguments
integer dimensions
integer shape(*)
double precision geometry(*)
integer i
do i=1,dimensions
GeomMin(i)=geometry(2*(i-1)+1)
GeomMax(i)=geometry(2*i)
enddo
return
end
c =====================================================
subroutine getbs(bl,sl)
c =====================================================
c returns the locations for the bubble and spike as
c computed from the initial perturbation
implicit none
c
include 'cuser.i'
double precision width
double precision bl,sl
c # width of the shock tube
width = GeomMax(2)-GeomMin(2)
c # used to convert modes to freq
bl = (width/omegaL)/4.d0
sl = (width/omegaL)
return
end
c =====================================================
subroutine spikeandbubble(np,scalar,xc,location)
c =====================================================
c given a vertical slice of data (a ray in the x direction)
c ordered from x-large to x-small find the first
c jump in the scalar ( from 0 to .95)
implicit none
integer np
double precision scalar(np)
double precision xc(np)
double precision location
integer i,ifirst,l
l=1
ifirst = 0
do i=1,np
if(scalar(i).gt.0.95) then
l=i
ifirst = 1
end if
if (ifirst.eq.1) go to 10
end do
10 location = xc(l)
return
end