Intensity-based Emission with pure absorption

Intensity-based (ibased) emission calculation with pure absorption is a method to compute the outgoing flux from the top of the atmosphere, with a given temperature-pressure and opacity profile, by transfering intensity (not flux) through the layers assuming no scattering.

\(I_0 (\mu) = I (\tau_B^\prime, \mu) e^{-\tau_B^\prime} +\int^{\tau_B^\prime}_{0} B (\tau) e^{-\tau^\prime} d \tau^\prime \approx B(T_B) e^{- \tau_{B}/\mu} + \sum_{n=0}^{N-1} B(T_n) (e^{-\tau_{n}/\mu} - e^{-\tau_{n+1}/\mu}).\)

For the \(\mathsf{N}\) -stream, the outgoing flux can be comuted by discretizing the flux using the Gauss-Legendre quadrature as

\(F_\mathrm{out} = 2 \pi \int_0^1 \mu I_0(\mu) d \mu \approx 2 \pi \sum_{i=1}^{\mathsf{N}^\prime} w_i \mu_i I_0 (\mu_i),\)

where \(\mathsf{N}^\prime = \mathsf{N}/2\) .

See Section 4.2 in Paper II for further details on the intensity-based emission calculation method.

Uses ArtPureEmis class

To perform tue ibased calculation using ArtEmisPure in ExoJAX, one simply needs to change the rtsolver option to either ibased or ibased_linsap.

The difference between ibased and ibased_linsap (ibased w/ linear source approximation) in ExoJAX lies in the assumptions regarding the distribution of the source function within the layer. The former assumes a uniform source function, while the latter assumes a linearly approximated source function. The latter can also be described as a linear version of Olson and Kunasz’s method.

In the case of ibased, the number of streams ( \(\mathsf{N}\) ) can be specified. This is done using the nstream option (which needs to be an even number).

from exojax.spec.atmrt import ArtEmisPure

art = ArtEmisPure(pressure_top=1.e-8,
                  pressure_btm=1.e2,
                  nlayer=100,
                  nu_grid=nu_grid,
                  rtsolver="ibased", # "ibased" or "ibased_linsap"
                  nstream=8)         # specify the number of the streams

art.change_temperature_range(400.0, 1500.0) #sets temperature range
Tarr = art.powerlaw_temperature(1300.0, 0.1) # sets a powerlaw T-P profile
mmr_arr = art.constant_mmr_profile(0.1) # sets a constant mass mixing ratio
gravity = art.constant_gravity_profile(2478.57) #sets a constant gravity

# we here call OpaPremodit as opa just to compute xsmatrix
opa = OpaPremodit(mdb=mdb,
                  nu_grid=nu_grid,
                  diffmode=diffmode,
                  auto_trange=[art.Tlow, art.Thigh])
xsmatrix = opa.xsmatrix(Tarr, art.pressure)

dtau = art.opacity_profile_xs(xsmatrix, mmr_arr, opa.mdb.molmass, gravity) # computes optical depth profile
F0 = art.run(dtau, Tarr) # computes spectrum

Layer-wise computation OpartEmisPure

You can further reduce device memory usage with layer-wise computation. Please use OpartEmisPure for this purpose.

Here is an example. The user needs to define a class that returns the optical depth profile for a given layer.

from exojax.spec.opacalc import OpaPremodit
from exojax.spec.opart import OpartEmisPure
from exojax.spec.layeropacity import single_layer_optical_depth
from exojax.utils.grids import wavenumber_grid
from exojax.spec.api import MdbExomol
from exojax.utils.astrofunc import gravity_jupiter
import jax.numpy as jnp
from jax import config
config.update("jax_enable_x64", True)

class OpaLayer:
    # user defined class, needs to define self.nugrid
    def __init__(self, Nnus=100000):
        self.nu_grid, self.wav, self.resolution = wavenumber_grid(
            #1900.0, 2300.0, Nnus, unit="cm-1", xsmode="premodit"
            2050.0, 2150.0, Nnus, unit="cm-1", xsmode="premodit"

        )
        self.mdb_co = MdbExomol(".database/CO/12C-16O/Li2015", nurange=self.nu_grid)
        self.opa_co = OpaPremodit(
            self.mdb_co,
            self.nu_grid,
            auto_trange=[500.0, 1500.0],
            dit_grid_resolution=1.0,
        )
        self.gravity = gravity_jupiter(1.0, 10.0)


    def __call__(self, params):
        temperature, pressure, dP, mixing_ratio = params
        xsv_co = self.opa_co.xsvector(temperature, pressure)
        dtau_co = single_layer_optical_depth(
            dP, xsv_co, mixing_ratio, self.mdb_co.molmass, self.gravity
        )
        return dtau_co

Then, the user can define the layer-wise computation as follows. However, users must define the layer’s update function themselves. This requirement is designed to prevent XLA compilation overhead.

opalayer = OpaLayer(Nnus=100000)
opart = OpartEmisPure(opalayer, pressure_top=1.0e-5, pressure_btm=1.0e1, nlayer=200, nstream=8)

def layer_update_function(carry_tauflux, params):
    carry_tauflux = opart.update_layer(carry_tauflux, params)
    return carry_tauflux, None

temperature = opart.clip_temperature(opart.powerlaw_temperature(1300.0, 0.1))
mixing_ratio = opart.constant_mmr_profile(0.01)
layer_params = [temperature, opart.pressure, opart.dParr, mixing_ratio]
flux = opart(layer_params, layer_update_function)

Refer to the Getting Started with Opart; GPU memory-efficient Emission Spectrum for another example of OpartEmisPure.