c =====================================================
subroutine combl()
c =====================================================
c
c Create and initialize application specific common-blocks.
c
implicit none
include 'cles.i'
include 'cuser.i'
include 'heairacesf6.i'
integer lin
double precision gamma1_hemx,gamma1_air,gamma1_sf6
double precision MachNo,gamma_hemx,rgas_hemx
double precision gm1, gp1, g2, M2
double precision c_Sair, c_Usf6
double precision c_hemx,c_air
double precision pwr
double precision vx(7), vec(7)
integer m, i
double precision x, th
integer nv
parameter( lin = 5 )
c
nv = 7
call stdcdat()
c
open(unit=lin,status='old',form='formatted',file='init.dat')
read (lin,* ) MachNo,dloc
read (lin, *) t_air,p_air
read (lin, *) t_hemx,fraction_air
read (lin, *) crtnL,crtnR
read (lin, *) zicki
c shock capturing info
c read (lin, *) dwidth, minc, minjump, gradS
read (lin, *) CurvP,CurvRho,CurvS
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 # 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 # speed of sounds
c_hemx = sqrt(gamma_hemx*rgas_hemx*t_hemx)
c_air = sqrt(gamma_air*rgas_air*t_air)
c # initially at rest
u_hemx = 0.d0
u_air = 0.d0
c # compute the pressure needed in He to make the required shock.
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)
rho_air = p_air/(rgas_air*t_air)
c # assume gas curtain (SF6) is at the same temp/pressure as air
p_sf6 = p_air
t_sf6 = t_air
rho_sf6 = p_sf6/(rgas_sf6*t_sf6)
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 # data in right state: -AIR
gamma1_air = gamma_air-1.d0
rhou_air = rho_air * u_air
e_air = p_air/gamma1_air + 0.5d0*rhou_air**2/rho_air
c
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 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)
endif
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