Viscous Heating Energy Conservation Issues

[astertaylor]

I've been trying to implement viscous heating in a very narrow cell region in a symmetric 2d spherical coordinate system. I've set the cell heating as I want, and I've calculated the total energy that should be provided by the heating. However, when I calculate the SED and integrate the luminosity in frequency, I find that the energy in the SED is ~1.8x the energy input by the cell heating. I've set the star temperature to be negligible, and the temperature distribution of the system is what I'd expect from the heating. I suspect there's some issue in the calculation of the SEDs from a system with heating. See a MWE for the problem_setup.py file below.

#
# Import NumPy for array handling
#
import numpy as np

def dust_density(r,theta,dtheta):
    if (theta+dtheta/2 >= np.pi/2)&(RX<=r<=RC):
        kappaP = kappa0 * kB*TX/(h*nu0) * 3.83223 * (r/RX)**(-3/4)
        return 1/(kappaP * r * np.sin(theta) * dtheta)
    else:
        return 0
    
def disk_heating(r,theta, dr, dtheta):
    if theta+dtheta <= np.pi/2:
        return 0
    elif (r > RC) or (r < RX): return 0
    else:
        return( σ*TX**4 * RX**3 /(r**4 * np.sin(theta) * dtheta) )
    
# Some natural constants
#
au  = 1.49598e13     # Astronomical Unit       [cm]
Mp  = 2e30     # Jupiter mass              [g]
Mstar = 2e33   # Solar mass              [g]
Rp  = 1e10        # Solar radius            [cm]

G = 6.67259e-8  # Gravitational constant   [cm3 g-1 s-2]
σ = 5.67051e-5  # Stefan-Boltzmann constant [erg cm-2 s-1 K-4]
c = 2.99792458e10 # Speed of light [cm/s]
h = 6.62607015e-27 # Planck constant [erg s]
kB = 1.380649e-16 # Boltzmann constant [erg K-1]

a = 5*au
Mpdot = 3* 6.02e16 # g / s

RH = a*(Mp/(3*Mstar))**(1/3)

nu0 = 1e14 
kappa0 = 10

#
# Monte Carlo parameters
#
nphot    = 1000000
#
# Grid parameters
#
nx       = 1000
ny       = 150        # One-sided only. So the "real" value is twice this.
nz       = 1
#
# Model parameters
#
rin      = Rp
rout     = RH
#
# Planet parameters
#
RC = RH/3
RX = 3.8e10

TX = 950

#
# Make the coordinates
#
xi       = rin * (rout/rin)**(np.linspace(0.,nx,nx+1)/(nx-1.0))
yi       = np.pi/2.0 * np.linspace(0.,1.,ny+1)
zi       = np.array([0.,2*np.pi])
xc       = 0.5e0 * ( xi[:-1] + xi[1:] )
yc       = 0.5e0 * ( yi[:-1] + yi[1:] )
#
# Make the dust density model
#
rr,tt    = np.meshgrid(xc,yc,indexing='ij')

rhod = np.zeros((nx, ny))
for i in range(nx):
    for j in range(ny):
        dtheta = yi[j+1] - yi[j]
        rhod[i,j] = dust_density(rr[i,j],tt[i,j], dtheta)
#
# Make the heating model
#
htot = 0
hdisk = np.zeros((nx, ny))
for i in range(nx):
    for j in range(ny):
        dr = xi[i+1] - xi[i]
        dtheta = yi[j+1] - yi[j]
        hdisk[i,j] = disk_heating(rr[i,j], tt[i,j], dr, dtheta)
        htot += hdisk[i,j] * rr[i,j]**2 * np.sin(tt[i,j]) * dr * dtheta * 2*np.pi
print(f"Total heating: {htot:0.2e} erg/s")
#
# Write the wavelength_micron.inp file
#
lam1     = 0.1e0
lam2     = 7.0e0
lam3     = 25.e0
lam4     = 1.0e4
n12      = 20
n23      = 100
n34      = 30
lam12    = np.logspace(np.log10(lam1),np.log10(lam2),n12,endpoint=False)
lam23    = np.logspace(np.log10(lam2),np.log10(lam3),n23,endpoint=False)
lam34    = np.logspace(np.log10(lam3),np.log10(lam4),n34,endpoint=True)
lam      = np.concatenate([lam12,lam23,lam34])
nlam     = lam.size
#
# Write the wavelength file
#
with open('wavelength_micron.inp','w+') as f:
    f.write('%d\n'%(nlam))
    for value in lam:
        f.write('%13.6e\n'%(value))
#
#
# Write the stars.inp file
#
with open('stars.inp','w+') as f:
    f.write('2\n')
    f.write('0 %d\n\n'%(nlam))
    # f.write('%13.6e %13.6e %13.6e %13.6e %13.6e\n\n'%(Rp*0.999999,Mp,pstar[0],pstar[1],pstar[2]))
    # for value in lam:
    #     f.write('%13.6e\n'%(value))
    # f.write('\n%13.6e\n'%(-Tp))
#
# Write the grid file
#
with open('amr_grid.inp','w+') as f:
    f.write('1\n')                       # iformat
    f.write('0\n')                       # AMR grid style  (0=regular grid, no AMR)
    f.write('100\n')                     # Coordinate system
    f.write('0\n')                       # gridinfo
    f.write('1 1 0\n')                   # Include x,y,z coordinate
    f.write('%d %d %d\n'%(nx,ny,nz))     # Size of grid
    for value in xi:
        f.write('%13.6e\n'%(value))      # X coordinates (cell walls)
    for value in yi:
        f.write('%13.6e\n'%(value))      # Y coordinates (cell walls)
    for value in zi:
        f.write('%13.6e\n'%(value))      # Z coordinates (cell walls)
#
# Write the density file
#
with open('dust_density.inp','w+') as f:
    f.write('1\n')                       # Format number
    f.write('%d\n'%(nx*ny*nz))           # Nr of cells
    f.write('1\n')                       # Nr of dust species
    data = rhod.ravel(order='F')         # Create a 1-D view, fortran-style indexing
    data.tofile(f, sep='\n', format="%13.6e")
    f.write('\n')
#
# Write the disk heating file
#
with open('heatsource.inp','w+') as f:
    f.write('1\n')                       # Format number
    f.write('%d\n'%(nx*ny*nz))           # Nr of cells
    data = hdisk.ravel(order='F')        # Create a 1-D view, fortran-style indexing
    data.tofile(f, sep='\n', format="%13.6e")
    f.write('\n')
#
# Dust opacity control file
#
with open('dustopac.inp','w+') as f:
    f.write('2               Format number of this file\n')
    f.write('1               Nr of dust species\n')
    f.write('============================================================================\n')
    f.write('1               Way in which this dust species is read\n')
    f.write('0               0=Thermal grain\n')
    f.write('linear        Extension of name of dustkappa_***.inp file\n')
    f.write('----------------------------------------------------------------------------\n')
#
# Write the radmc3d.inp control file
#
with open('radmc3d.inp','w+') as f:
    f.write('nphot = %d\n'%(nphot))
    f.write('scattering_mode_max = 0\n')   # Put this to 1 for isotropic scattering

#
# Write the linear opacity file
#
with open('dustkappa_linear.inp','w+') as f:
    f.write('1\n')                       # Format number
    f.write('%d\n'%(1000))
    for value in np.logspace(-1, 4, 1000):
        f.write('%13.6e  %13.6e\n'%(value, kappa0*1e4*(c/(nu0*value))))