c
c Copyright (C) 2002 Ralf Deiterding
c Brandenburgische Universitaet Cottbus
c
c Copyright (C) 2003-2007 California Institute of Technology
c Ralf Deiterding, ralf@amroc.net
c
c # Calculate detonation profile by integrating
c # production rates in space. Detonation is assumed
c # to be traveling from left to right.
c =====================================================
subroutine ic(maxmx,meqn,mbc,mx,x,dx,q)
c =====================================================
implicit double precision (a-h,o-z)
c
common /grkstat/ jsuc,jfail,jf,ja
include "ck.i"
include "cuser.i"
c
dimension q(meqn,1-mbc:maxmx+mbc)
dimension x(1-mbc:maxmx+mbc)
parameter ( lengrkaux = 1 )
dimension grkaux(lengrkaux), scale(LeNsp), qint(LeNsp)
external fcn, Jfcn, odeout
c
tol = 1.d-05
sclimit = 1.d-03
uround = 5.d-15
tst = dx
ts = 0.d0
te = 0.d0
tout = 0.d0
ichan = 0
ijac = 0
iscale = 0
jsuc = 0
jfail = 0
jf = 0
ja = 0
Nint = Nsp
hmax = dx
c
do k=1, Nsp
qint(k) = Yk(k)
enddo
c
c # Integrate production function up to state necessary
c # on right side of current grid.
if (x(mx) .lt. sloc) then
te = sloc-x(mx)
call grk4a(Nint,fcn,Jfcn,ts,qint,te,tol,
& hmax,tst,lengrkaux,grkaux,ichan,ijac,odeout,
& iscale,scale,lengrkw,grkwork,lenint,igrkwork)
endif
c
c # Fill current grid from right to left and integrate
c # production function by going from cell to cell.
do 150 i=mx+mbc,1-mbc,-1
if (x(i) .ge. sloc) then
do k=1, Nsp
q(k,i) = rho0*Yk(k)
enddo
q(Nsp+1,i) = rhou0
q(Nsp+2,i) = E0
q(Nsp+3,i) = T0
else
ts = te
te = sloc-x(i)
call grk4a(Nint,fcn,Jfcn,ts,qint,te,tol,
& hmax,tst,lengrkaux,grkaux,ichan,ijac,odeout,
& iscale,scale,lengrkw,grkwork,lenint,igrkwork)
call RankineH(Nsp,qint,udet,rho,u,p,T,ifail)
do k=1, Nsp
q(k,i) = rho*qint(k)
enddo
u = udet - u
if (Nstat.eq.1) u = u - udet
q(Nsp+1,i) = rho*u
q(Nsp+3,i) = T
q(Nsp+2,i) = rho * (0.5d0*u**2 +
& avgtabip(T,qint,hms,Nsp)) - p
endif
q(Nsp+4,i) = 0.d0
c
150 continue
c
do 160 i=1,mx
if (x(i).lt.x1p.or.x(i).gt.x2p) goto 160
rhoo = 0.d0
rho = 0.d0
rhoW = 0.d0
do k=1,Nsp
rhoo = rhoo + q(k,i)
q(k,i) = rhop*Ykp(k)
rho = rho + q(k,i)
rhoW = rhoW + q(k,i)/Wk(k)
enddo
q(Nsp+3,i) = Tp
u = q(Nsp+1,i)/rhoo
q(Nsp+1,i) = u*rho
p = rhoW*RU*q(Nsp+3,i)
q(Nsp+2,i) = 0.5d0*(q(Nsp+1,i)**2)/rho
& + avgtabip(q(Nsp+3,i),q(1,i),hms,Nsp)-p
160 continue
c
return
end
c
c # For given mass fraction evaluate data on
c # Hugoniot curve and return production rates in
c # mass fraction. Note that this integration is in
c # space.
c =========================================================
subroutine fcn(n,tt,Y,dY,na,aux,ifcnfail)
c =========================================================
implicit double precision (a-h,o-z)
include "ck.i"
include "cuser.i"
dimension Y(n),dY(n),aux(na),Ck(LeNsp)
c
ifcnfail = 0
call RankineH(n,Y,udet,rho,u,p,T,ifcnfail)
do k = 1, Nsp
Ck(k) = rho*Y(k)/Wk(k)
enddo
call prodckwc(T,Ck,dY,n)
do k = 1,Nsp
dY(k) = dY(k)*Wk(k)/(rho*u)
enddo
c
return
end
c
c # The Jacobian is approximated and therefore left empty.
c =========================================================
subroutine Jfcn(n,tt,df,ddf,na,aux,ifcnfail)
c =========================================================
implicit double precision (a-h,o-z)
dimension df(n),ddf(n,n),aux(na)
c
ifcnfail = 0
c
return
end
c
c =========================================================
subroutine odeout(n,tt,f,naux,aux,ind)
c =========================================================
implicit double precision(a-h,o-z)
include "ck.i"
dimension f(n), aux(naux)
common /grkstat/ jsuc,jfail,jf,ja
common /output/ tout
c
write (6,*) tt+tout,jfail,f(1)
return
end