Note

Download full notebook here

Adding a Variable \(f_\text{sed}\)

In this tutorial you will:

1. Learn how to use variable \(f_\text{sed}\) in your calculation

2. Compare outputs of constant and variable \(f_\text{sed}\)

Citation for this methodology: https://ui.adsabs.harvard.edu/abs/2021arXiv211005903R/abstract

import warnings
warnings.filterwarnings('ignore')

from bokeh.io import output_notebook
from bokeh.plotting import show, figure
from bokeh.palettes import Colorblind
output_notebook()
import numpy as np
import pandas as pd
import astropy.units as u
import time

#here is pyeddy
import virga.justdoit as jdi
import virga.justplotit as jpi

Loading BokehJS ...

Let us again load in a Hot Jupiter PT profile and run virga for constant \(f_\text{sed}=1\)

mieff_directory = "/data/virga/"
metallicity = 1 #atmospheric metallicity relative to Solar
mean_molecular_weight = 2.2 # atmospheric mean molecular weight

#set the run
a = jdi.Atmosphere(['MnS','Cr','MgSiO3','Fe'],
                  fsed=1,mh=metallicity,
                 mmw = mean_molecular_weight)

#set the planet gravity
a.gravity(gravity=7.460, gravity_unit=u.Unit('m/(s**2)'))

#Get preset pt profile for testing
a.ptk(df = jdi.hot_jupiter())

#get full dictionary output
all_out = jdi.compute(a, as_dict=True, directory=mieff_directory)

Not doing sublayer as cloud deck at the bottom of pressure grid

Variable \(f_\text{sed}\) expression

We have derived an altitude-dependent expression for fsed, given by

\(f_\text{sed}(z) = \alpha\exp\left(\frac{z-z_*}{6\beta H_0}\right) + \epsilon\)

where \(\alpha\) and \(\beta\) are free parameters, \(z_*\) is the altitude at which \(f_\text{sed}(z_*)=\alpha+\epsilon\) (default=top of atmosphere), \(H_0\) is the scale-height at the effective temperature of the planet and \(\epsilon\) is the minimum value of \(f_\text{sed}\) we wish to consider (default=1e-2).

Varying \(\alpha\) controls the value of \(f_\text{sed}\) at \(z_*\) (here the top of the atmosphere)

pressure = all_out['pressure']
z = all_out['altitude']
H = all_out['scale_height']
alpha = [1,10,20,100]; beta = [0.5]
fig = jpi.plot_fsed(pressure, z, H, alpha, beta)
fig.legend.location='top_left'
show(fig)

Varying \(\beta\) controls the slope of \(f_\text{sed}\), or the rate of change of \(f_\text{sed}\).

alpha = [10]; beta = [0.25,0.5,1,5]
fig = jpi.plot_fsed(pressure, z, H, alpha, beta)
fig.legend.location='top_left'
show(fig)

## Running code with variable \(f_\text{sed}\)

To run VIRGA with this variable \(f_\text{sed}\) expression, we must set the parameterization to “exp” and define fsed $ =:nbsphinx-math:alpha`$ and b $ =:nbsphinx-math:beta`$ in the Atmosphere class.

# define alpha and beta
alpha = 10 # generally choose alpha to be between 1 and 100
beta = 0.5 # generally choose beta to be between 0.1 and 5

#set the run
a = jdi.Atmosphere(['MnS','Cr','MgSiO3','Fe'],
                   mh=metallicity, mmw = mean_molecular_weight,
                   fsed=alpha, b=beta, param="exp")

#set the planet gravity
a.gravity(gravity=7.460, gravity_unit=u.Unit('m/(s**2)'))

#Get preset pt profile for testing
a.ptk(df = jdi.hot_jupiter())

#get full dictionary output
all_out_var = jdi.compute(a, as_dict=True, directory=mieff_directory)

Not doing sublayer as cloud deck at the bottom of pressure grid

We can also visualise the fsed used in the virga calculation using jpi.fsed_from_output

fig = jpi.fsed_from_output(all_out_var,['variable fsed'])
fig.legend.visible=False
show(fig)

Additional flexibility

You can also change the default parameters \(z_*\), \(\epsilon\) and \(H_0\).

alpha=10;
beta=0.5;
eps=1e-1; # this stops fsed being smaller than 1e-1
alpha_pressure=1e-4; # z_* will be set to the altitude at 1e-4 bars
Teff=2000; # scale-height H_0 will be defined at this temperature

# change alpha, beta, epsilon and parameterisation
a = jdi.Atmosphere(['MnS','Cr','MgSiO3','Fe'],
                   mh=metallicity, mmw = mean_molecular_weight,
                   fsed=alpha, b=beta, param='exp',
                   eps=eps)

a.gravity(gravity=7.460, gravity_unit=u.Unit('m/(s**2)'))

# change alpha_pressure and Teff
a.ptk(df = jdi.hot_jupiter(), alpha_pressure=alpha_pressure, Teff=Teff)

Compare constant to variable \(f_\text{sed}\) calculations

#   compare solvers
labels = ["constant", "variable"]
param = ['const', 'exp']
fsed = [1, 10] # fsed=1 for constant, alpha=10 for variable
beta = 0.5 # note b is irrelevant when param='const'

output = []
for i, p, fs in zip(range(2), param, fsed):
    print(labels[i] + " fsed")
    t1 = time.time()
    a = jdi.Atmosphere(['MnS','Cr','MgSiO3','Fe'],
                   mh=metallicity, mmw = mean_molecular_weight,
                   fsed=fs, b=beta, param=p)

    a.gravity(gravity=7.460, gravity_unit=u.Unit('m/(s**2)'))
    a.ptk(df = jdi.hot_jupiter())

    all_out = jdi.compute(a, as_dict=True, directory=mieff_directory)
    t2 = time.time() - t1
    print("time = %.2f seconds" % t2)
    output.append(all_out)

constant fsed
Not doing sublayer as cloud deck at the bottom of pressure grid
time = 3.30 seconds
variable fsed
Not doing sublayer as cloud deck at the bottom of pressure grid
time = 3.78 seconds

Introducing an altitude-dependent expression for \(f_\text{sed}\) does not have detrimental effects on the computational efficiency of VIRGA.

We can plot the mean mass mixing ratio of the condensates for each solver. We use a solid line for constant \(f_\text{sed}\) and a dashed line for the variable \(f_\text{sed}\).

fig = jpi.condensate_mmr(output, plot_height=400, plot_width=500,x_range=[1e-9,1e-2])
fig.legend.location='bottom_left'
show(fig)

The comparison plot shows the results for all of the gases for both solvers. For a clearer comparison, you can choose which gas to view.

show(jpi.condensate_mmr(output, gas='Fe'))

We can similarly analyse the particle size distributions produced by the two solvers. Here we plot 1. the mean particle radius for each condensate with pressure 2. the particle distribution at a given pressure level (1e-3)

fig, dndr = jpi.radii(output,at_pressure=1e-3)
show(fig)

Again, you can choose which gas to view.

fig, dndr = jpi.radii(output,gas="MnS",at_pressure=1e-3)
show(fig)

show(jpi.opd_by_gas(output))
show(jpi.opd_by_gas(output, gas="MgSiO3"))

Finally we can compare the cloud optical depth, single scattering and asymmetry parameters.

for i in range(len(output)):
    print(labels[i] + " solver")
    show(jpi.all_optics(output[i]))

constant solver

variable solver