[FLASH-BUGS] small bug in conductivity
David J. Lin
dlin at vela.astro.northwestern.edu
Tue Apr 2 13:17:24 CST 2002
Hi,
In the current conductivity.F90 for the "stellar" conductivity module
(/source/materials/conductivity/stellar/), there appears to be a bug which
is related to the recent switch from FLASH1.6 to FLASH2.0. Specifically,
the current routine wants to include "network_common.fh", which is no
longer used in FLASH2.0.
I've attached the file with my simple fixes to the problem, calling the
multifluid database to fill in the zion and aion variables. For your
convenience, I prefaced with "### NU ADD ###" the lines which I added or
changed.
I'm fairly certain these fixes will do the trick, but if I've made a
mistake, or if there's an easier or more efficient fix, please let me
know.
Thanks.
David
-------------- next part --------------
subroutine conductivity(xtemp,xden,massfrac,cond,diff_coeff)
use dBase, ONLY: ionmax
use ModuleEos, ONLY: eos
! ### NU ADD #### 04.02.02 DJL
use multifluid_database, ONLY: get_mfluid_property
! ### ### ###
implicit none
save
!
! this is Frank's routine that approximates the thermal transport
! coefficients for electron conduction.
!
! The calling sequence has been modified for use in FLASH. Instead of
! passing in the electron pressure and degeneracy pressure, this routine
! makes an eos call to the general FLASH eos driver (eos_fcn). This
! makes use of abar and zbar which were computed in this routine anyway.
! The mass fractions are passed in through the array massfrac instead of
! xmass, to avoid using the common blocks. The only information which
! comes from the common blocks is the list of isotope proton numbers and
! atomic masses (aion and zion) and the number of isotopes.
!
! input:
! xtemp = temperature temp (in K)
! xden = density den (in g/cm**3)
! massfrac = mass fractions of the composition
!
! output:
! cond = conductivity
! diff_coeff = diffusion coefficient ( = cond/(rho*c_v))
!
!
! modified 7-9-99
!
! now returns the conductivity, not opacity
! eliminated zion, aion, and ionmax from the argument list,
! these come from the network includes
!
! ### NU ADD ### 4.2.02. DJL
! Commented out the following ... seems to be a throwback from FLASH1.6
! include 'network_common.fh'
! ### ### ###
! ### NU ADD ### 4.2.02. DJL
! Add declaration for zion and aion
real, save :: zion(ionmax), aion(ionmax)
! ### ### ###
! declare the pass
real xtemp, xden, massfrac(ionmax), cond, diff_coeff
! declare the variables used by the eos
real pres, ener, dpt, dpd, det, ded, gammac
! declare the electron contributions -- these are used in this routine
real pep, xne, eta
real opac
real c_p, chid
!..declare the internal variables
integer iz,i
real xmu,t6,orad,ocond,ytot1,ymassfrac,abar,zbar,w(6),xh, &
& xhe,xz,xkc,xkap,xkb,xka,dbar,oiben1,d0log,xka1, &
& xkw,xkaz,dbar1log,dbar2log,oiben2,t4,t4r,t44,t45,t46, &
& ck1,ck3,ck2,ck4,ck5,ck6,xkcx,xkcy,xkcz,ochrs,th,fact, &
& facetax,faceta,ocompt,tcut,cutfac,xkf,dlog10,zdel,zdell10, &
& eta0,eta02,thpl,thpla,cfac1,cfac2,oh,pefac,pefacl,pefacal, &
& dnefac,wpar2,walf,walf10,thx,thy,thc,farg,ffac,xmas,ymas, &
& wfac,cint,vie,cie,tpe,yg,xrel,beta2,jy,vee,cee,ov1,ov, &
& drel,drel10,drelim
! various physical and derived constants
! con2 = con1*sqrt(4*pi*e*e/me)
! meff = hbar/(me*c)*(3*pi**2)**(1/3)
! weid = (pi*kerg)**2/(3*me)
! iec = 4*e**4*me/(3*pi*hbar**3)
! xec = hbar/kerg*e*sqrt(4*pi/me)
!
real third,twoth,pi,avo,c,ssol,asol,zbound,t7peek, &
& con2,k2c,meff,weid,iec,xec,rt3
parameter (third = 1.0e0/3.0e0, &
& twoth = 2.0e0 * third, &
& pi = 3.1415926535897932384e0, &
& avo = 6.0221367e23, &
& c = 2.99792458e10, &
& ssol = 5.67050407222e-5, &
& asol = 4.0e0*ssol/c, &
& zbound = 0.1e0, &
& t7peek = 1.0e20, &
& k2c = 4.0e0/3.0e0*asol*c, &
& meff = 1.194648642401440e-10, &
& weid = 6.884326138694269e-5, &
& iec = 1.754582332329132e16, &
& xec = 4.309054377592449e-7, &
& rt3 = 1.7320508075688772e0, &
& con2 = 1.07726359439811217e-7)
! conversion coefficient for opacity to conductivity
real cond_factor
parameter (cond_factor = 4.e0*asol*c/3.e0)
! ### NU ADD ### 04.02.02 DJL
! First call - fill in zion and aion
logical, save :: firstCall = .TRUE.
!===========================================================================
if (firstCall) then
! Define zion and aion via multifluid database
call get_mfluid_property ("num positive", zion)
call get_mfluid_property ("num total", aion)
firstCall = .FALSE.
endif
!===========================================================================
! ### ### ###
!..initialize
opac = 0.0e0
orad = 0.0e0
ocond = 0.0e0
oiben1 = 0.0e0
oiben2 = 0.0e0
ochrs = 0.0e0
oh = 0.0e0
ov = 0.0e0
zbar = 0.0e0
ytot1 = 0.0e0
!..set the composition variables
do i=1,6
w(i) = 0.0e0
enddo
do i = 1,ionmax
iz = min(3,max(1,int(zion(i))))
ymassfrac = massfrac(i)/aion(i)
w(iz) = w(iz) + massfrac(i)
w(iz+3) = w(iz+3) + zion(i) * zion(i) * ymassfrac
zbar = zbar + zion(i) * ymassfrac
ytot1 = ytot1 + ymassfrac
enddo
abar = 1.0e0/ytot1
zbar = zbar * abar
t6 = xtemp * 1.0e-6
xh = w(1)
xhe = w(2)
xz = w(3)
! now that we have abar and zbar, call the eos and get the electron
! contributions
call eos(xden,xtemp,pres,ener,massfrac,abar,zbar, &
& dpt,dpd,det,ded,gammac,pep,xne,eta,1)
! now generate the specific heat at constant pressure from these quantities
! see C&G 9.98
chid = dpd*xden/pres
c_p = det*gammac/chid
!..radiative section:
!..from iben apj 196 525 1975
if (xh .lt. 1.0e-5) then
xmu = max(1.0e-99, w(4)+w(5)+w(6)-1.0e0)
xkc = (2.019e-4*xden/t6**1.7e0)**(2.425e0)
xkap = 1.0e0 + xkc * (1.0e0 + xkc/24.55e0)
xkb = 3.86e0 + 0.252e0*sqrt(xmu) + 0.018e0*xmu
xka = 3.437e0 * (1.25e0 + 0.488e0*sqrt(xmu) + 0.092e0*xmu)
dbar = exp(-xka + xkb*log(t6))
oiben1 = xkap * (xden/dbar)**(0.67e0)
end if
if ( .not.((xh.ge.1.0e-5) .and. (t6.lt.1.0)) .and. &
& .not.((xh.lt.1.0e-5) .and. (xz.gt.zbound)) ) then
if (t6 .gt. 1.0) then
d0log = -(3.868e0 + 0.806e0*xh) + 1.8e0*log(t6)
else
d0log = -(3.868e0 + 0.806e0*xh) + &
& (3.42e0 - 0.52e0*xh)*log(t6)
endif
xka1 = 2.809e0 * exp(-(1.74e0 - 0.755e0*xh) &
& * (log10(t6) - 0.22e0 + 0.1375e0*xh)**2)
xkw = 4.05e0 * exp(-(0.306e0 - 0.04125e0*xh) &
& * (log10(t6) - 0.18e0 + 0.1625e0*xh)**2)
xkaz = 50.0e0*xz*xka1 * exp(-0.5206e0* &
& ((log(xden)-d0log)/xkw)**2)
dbar2log = -(4.283e0 + 0.7196e0*xh) + 3.86e0*log(t6)
dbar1log = -5.296e0 + 4.833e0*log(t6)
if (dbar2log .lt. dbar1log) dbar1log = dbar2log
oiben2 = (xden/exp(dbar1log))**(0.67e0) * exp(xkaz)
end if
!..from christy apj 144 108 1966
if ((t6.lt.1.5) .and. (xh .ge. 1.0e-5)) then
t4 = xtemp * 1.0e-4
t4r = sqrt(t4)
t44 = t4**4
t45 = t44 * t4
t46 = t45 * t4
ck1 = 2.0e6/t44 + 2.1e0*t46
ck3 = 4.0e-3/t44 + 2.0e-4/xden**(0.25e0)
ck2 = 4.5e0*t46 + 1.0e0/(t4*ck3)
ck4 = 1.4e3*t4 + t46
ck5 = 1.0e6 + 0.1e0*t46
ck6 = 20.0e0*t4 + 5.0e0*t44 + t45
xkcx = xh*(t4r/ck1 + 1.0e0/ck2)
xkcy = xhe*(1.0e0/ck4 + 1.5e0/ck5)
xkcz = xz*(t4r/ck6)
ochrs = pep * (xkcx + xkcy + xkcz)
end if
!..opacity in presence of hydrogen
if (xh .ge. 1.0e-5) then
if (t6 .lt. 1.0) then
orad = ochrs
else if (t6 .le. 1.5) then
orad = 2.0e0*(ochrs*(1.5e0 - t6) + oiben2*(t6 - 1.0e0))
else
orad = oiben2
end if
!..opacity in absence of hydrogen
else
if (xz .gt. zbound) then
orad = oiben1
else
orad = oiben1*(xz/zbound) + oiben2*((zbound-xz)/zbound)
end if
end if
!..add in the compton scattering opacity, weaver et al. apj 1978 225 1021
th = min(511.0e0, xtemp * 8.617e-8)
fact = 1.0e0 + 2.75e-2*th - 4.88e-5*th*th
facetax = 1.0e100
if (eta .le. 500.0) facetax = exp(0.522e0*eta - 1.563e0)
faceta = 1.0e0 + facetax
ocompt = 6.65205e-25/(fact * faceta) * xne/xden
orad = orad + ocompt
!..cutoff radiative opacity when 4kt/hbar is less than the plasma frequency
tcut = con2 * sqrt(xne)
if (xtemp .lt. tcut) then
if (tcut .gt. 200.0*xtemp) then
orad = orad * 2.658e86
else
cutfac = exp(tcut/xtemp - 1.0e0)
orad = orad * cutfac
end if
end if
!..fudge molecular opacity for low temps
xkf = t7peek * xden * (xtemp * 1.0e-7)**4
orad = xkf * orad/(xkf + orad)
!..conductivity section:
!..drel is the dividing line between nondegenerate and degenerate regions,
!..taken from clayton eq. 2-34. if the density is larger than drel, then
!..use the degenerate expressions. if the density is smaller than
!..drelim, use the non-degenerate formulas. in between drel and drelim,
!..apply a smooth blending of the two.
drel = 2.4e-8 * zbar/abar * xtemp * sqrt(xtemp)
if (xtemp .le. 1.0e5) drel = drel * 15.0e0
drel10 = log10(drel)
drelim = drel10 + 0.3e0
dlog10 = log10(xden)
!..from iben apj 196 525 1975 for non-degenerate regimes
if (dlog10 .lt. drelim) then
zdel = xne/(avo*t6*sqrt(t6))
zdell10 = log10(zdel)
eta0 = exp(-1.20322e0 + twoth * log(zdel))
eta02 = eta0*eta0
!..thpl factor
if (zdell10 .lt. 0.645) then
thpl = -7.5668e0 + log(zdel * (1.0e0 + 0.024417e0*zdel))
else
if (zdell10 .lt. 2.5) then
thpl = -7.58110e0 + log(zdel*(1.0e0 + 0.02804e0*zdel))
if (zdell10 .ge. 2.0) then
thpla = thpl
thpl = -11.0742e0 + &
& log(zdel**2 * (1.0e0 + 9.376e0/eta02))
thpl = 2.0e0*((2.5e0-zdell10)*thpla + &
& (zdell10-2.0e0)*thpl)
end if
else
thpl = -11.0742e0 + &
& log(zdel**2 * (1.0e0 + 9.376e0/eta02))
end if
end if
!..pefac and walf factors
if (zdell10 .lt. 2.0) then
pefac = 1.0e0 + 0.021876e0*zdel
if (zdell10 .gt. 1.5) then
pefacal = log(pefac)
pefacl = log(0.4e0 * eta0 + 1.64496e0/eta0)
cfac1 = 2.0e0 - zdell10
cfac2 = zdell10 - 1.5e0
pefac = exp(2.0e0 * (cfac1*pefacal + cfac2*pefacl))
end if
else
pefac = 0.4e0 * eta0 + 1.64496e0/eta0
end if
if (zdel.lt.40.0) then
dnefac = 1.0e0 + zdel * (3.4838e-4 * zdel - 2.8966e-2)
else
dnefac = 1.5e0/eta0 * (1.0e0 - 0.8225e0/eta02)
endif
wpar2 = 9.24735e-3 * zdel * &
& (xden*avo*(w(4)+w(5)+w(6))/xne + dnefac)/(sqrt(t6)*pefac)
walf = 0.5e0 * log(wpar2)
walf10 = 0.5e0 * log10(wpar2)
!..thx, thy and thc factors
if (walf10 .le. -3.0) then
thx = exp(2.413e0 - 0.124e0*walf)
else if (walf10 .le. -1.0) then
thx = exp(0.299e0 - walf*(0.745e0 + 0.0456e0*walf))
else
thx = exp(0.426e0 - 0.558e0*walf)
end if
if (walf10 .le. -3.0) then
thy = exp(2.158e0 - 0.111e0*walf)
else if (walf10 .le. 0.0) then
thy = exp(0.553e0 - walf*(0.55e0 + 0.0299e0*walf))
else
thy = exp(0.553e0 - 0.6e0*walf)
end if
if (walf10 .le. -2.5) then
thc = exp(2.924e0 - 0.1e0*walf)
else if (walf10 .le. 0.5) then
thc = exp(1.6740e0 - walf*(0.511e0 + 0.0338e0*walf))
else
thc = exp(1.941e0 - 0.785e0*walf)
end if
oh = (xh*thx + xhe*thy + w(6)*third*thc) / (t6*exp(thpl))
end if
!..from yakovlev & urpin soviet astro 1980 24 303 and
!..potekhin et al. 1997 aa 323 415 for degenerate regimes
if (dlog10 .gt. drel10) then
xmas = meff * xne**third
ymas = sqrt(1.0e0 + xmas*xmas)
wfac = weid * xtemp/ymas * xne
cint = 1.0e0
!..ion-electron collision frequency and the thermal conductivity
vie = iec * zbar * ymas * cint
cie = wfac/vie
!..electron-electron collision frequency and thermal conductivity
tpe = xec * sqrt(xne/ymas)
yg = rt3 * tpe/xtemp
xrel = 1.009e0 * (zbar/abar * xden * 1.0e-6)**third
beta2 = xrel**2/(1.0e0 + xrel**2)
jy = (1.0e0 + 6.0e0/(5.0e0*xrel**2) + 2.0e0/(5.0e0*xrel**4)) &
& * ( yg**3 / (3.0e0 * (1.0e0 + 0.07414 * yg)**3) &
& * log((2.81e0 - 0.810*beta2 + yg)/yg) &
& + pi**5/6.0e0 * (yg/(13.91e0 + yg))**4 )
vee = 0.511e0 * xtemp**2 * xmas/ymas**2 * sqrt(xmas/ymas) * jy
cee = wfac/vee
!..total electron thermal conductivity and conversion to an opacity
ov1 = cie * cee/(cee + cie)
ov = k2c/(ov1*xden) * xtemp**3
end if
!..blend the opacities in the intermediate region
if (dlog10 .le. drel10) then
ocond = oh
else if (dlog10 .gt. drel10 .and. dlog10 .lt. drelim) then
farg = pi * (dlog10 - drel10) / 0.3e0
ffac = 0.5e0 * (1.0e0 - cos(farg))
ocond = exp((1.0e0-ffac)*log(oh) + ffac*log(ov))
else if (dlog10 .ge. drelim) then
ocond = ov
end if
!..total opacity
opac = orad * ocond / (ocond + orad)
! now return the conductivity
cond = cond_factor*xtemp**3/(opac*xden)
diff_coeff = cond/(xden*c_p)
return
end
More information about the flash-bugs
mailing list