Note

Download full notebook here

Switching the Solver

In this tutorial you will learn:

1. How to switch solvers for calculating mass mixing ratios,

2. How to switch solvers for calculating particle size distributions.

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.

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(),latent_heat=True)

#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

The standard code uses an analytical solution to the eddy diffusion/gravitational settling equation for the mass mixing ratio (equation (4) in Ackerman & Marley 2001). However, an alternative approach is to solve this equation directly. This allows for features such as pressure-dependent \(f_{sed}\) and additional source terms.

We will now run the code for the original, analytical solution as well as the new, direct solver to verify the agreement between the results.

Note: directly solving the ODE is understandably slower than an analytical solution, therefore the new solver is more expensive than the original solver

#   compare solvers
labels = ["original", "direct"]
og_solver = [True, False]  # true means original, false is new solver
output = []
for i in range(2):
    print(labels[i] + " solver")
    t1 = time.time()
    all_out = jdi.compute(a, as_dict=True, directory=mieff_directory, og_solver=og_solver[i])
    t2 = time.time() - t1
    print("time = %.2f seconds" % t2)
    output.append(all_out)

original solver
Not doing sublayer as cloud deck at the bottom of pressure grid
time = 4.20 seconds
direct solver
Not doing sublayer as cloud deck at the bottom of pressure grid
time = 13.18 seconds

Mean Mass Mixing Ratio of the Condensates

We can plot the mean mass mixing ratio of the condensates for each solver. We use a solid line for the original, analytical solution and a dashed line for the new, direct solution.

show(jpi.condensate_mmr(output))

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'))

Particle Radii

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)

output[0].keys()

dict_keys(['pressure', 'pressure_unit', 'temperature', 'temperature_unit', 'wave', 'wave_unit', 'condensate_mmr', 'cond_plus_gas_mmr', 'mean_particle_r', 'droplet_eff_r', 'r_units', 'column_density', 'column_density_unit', 'opd_per_layer', 'single_scattering', 'asymmetry', 'opd_by_gas', 'condensibles', 'scalar_inputs', 'fsed', 'altitude', 'layer_thickness', 'z_unit', 'mixing_length', 'mixing_length_unit', 'kz', 'kz_unit', 'scale_height', 'cloud_deck'])

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)

Cumulative Optical Depth By Gas

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

Visualize Optical Depth, Asymmetry, and Single Scattering

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]))

original solver

direct solver

Calculating particle distributions

Currently, all particle distributions in virga are log-normal. The original solver utilises analytical expressions for the mean radius \(r_g\), effective radius \(r_{eff}\), and number concentration \(N\) given in Ackerman & Marley 2001.

However, if we were to extend the model using a different particle distribution that is not log-normal, these expressions would no longer be valid.

We have therefore introduced an alternative solver for direct computation of \(r_g\), \(r_{eff}\) and \(N\) which allows for general particle distributions.

t1 = time.time()
# to use the direct solver for particle distribution, set analytical_rg = False
all_out = jdi.compute(a, as_dict=True, directory=mieff_directory,
                      og_solver=False, analytical_rg=False)
output.append(all_out)
t2 = time.time() - t1
print("time = %.2f seconds" % t2)

Not doing sublayer as cloud deck at the bottom of pressure grid
time = 25.96 seconds

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

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

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

# all optics
labels.append("direct (+ direct size distn)")
for i in range(len(output)):
    print(labels[i] + " solver")
    show(jpi.all_optics(output[i]))

original solver

direct solver

direct (+ direct size distn) solver