Pyharp Documentation

Contents

• Working with the opacity class
  • Example 1. Compute Sonora2020 molecular opacities
  • Example 2. Compute Hydrogen continuum opacities
  • Example 3. Add a new opacity
  • Summary
• Built-in opacity sources
  • Sonora2020 molecular opacities
  • Hydrogen and Helium continuum
  • NetCDF dump-backed molecular opacities
  • References
• Working with the radiation class
• API Reference
  • Global States
  • Radiation Classes
  • Opacity Classes
  • Helper Functions
• pyharp-plot CLI
  • Basic Usage
  • Subcommands
  • Shared Options
  • Target Selection
  • Broadening Gases
  • Fallback Behavior
  • Output Examples
• pyharp-dump CLI
  • Basic Usage
  • Subcommands
  • Shared Options
  • Target Selection
  • Broadening Behavior
  • Multi-band Output
  • State Grid Output
  • NetCDF Naming Conventions
  • Global Attributes
  • Examples
• Troubleshooting

Introduction

Pyharp module provides a python interface to the C++ version of the HARP (High-performance Atmospheric Radiation Package) program [1]. It contains a minimum set of classes and subroutines for computing the radiative fluxes and radiances for a plane-parallel atmosphere. Pyharp is designed with efficiency – all heavy computations are implemented in C/C++ and wraped with a thin python interface. It can be used standalone to compute multi-column radiative transfer or as a module in other larger programs, both in C/C++ and Python. Significant changes have been made to the original HARP code as the backend data structure has changed from Athena to PyTorch.

Pyharp treats the radiative transfer problem as layers of computation, akin to a neural network, where the first layer – the input – is the atmospheric state; the second layer is the optical properties, and the final layer – the output – is the radiative fluxes or radiances. Each layer is represented by a class, which contains one important method: forward. The forward method takes torch.Tensor data as inputs and also returns torch.Tensor data as outputs. This design allows for chaining multiple layers together to perform a sophisticated computation as well as easy integration with other modules, such as neural networks or JIT (Just-In-Time) compilation.

Having both input and output as tensors makes the data flow through the layers transparent, mimicking the paradigm of functional programming.

Automatic differentiation is available when the input tensor has the torch.Tensor.requires_grad attribute set to True. This allows for the computation of gradients with respect to the input tensor, which is useful for optimization tasks.

The following example gives a practical demonstration of how to use the Pyharp module.

Example of calculating the radiative fluxes of Jupiter

Radiation configuration file

The program starts with writing a radiation configuration file in YAML format. Radiative transfer is a physical calculation regarding the interaction of radiation with molecule and matter. We specify the species, their composition and radiation bands in a configuration file that looks like this:

species:
  - name: H2mix
    composition: {H: 1.9, He: 0.1}

opacities:
  H2-molecule:
    type: multiband-ck
    data: ["sonora_2020_feh+000_co_100.data.196.pt"]
    species: [H2mix]

  H2-continuum:
    type: wavetemp
    fractions: [0.9, 0.1]
    data: ["H2-H2-eq.xiz.pt", "H2-He-eq.xiz.pt"]
    species: [H2mix]

bands: [sonora196]

sonora196:
  range: [30.8, 38300.] # wavenumber
  opacities: [H2-molecule, H2-continuum]
  solver: disort
  integration: weight
  flags: lamber,quiet,onlyfl,planck

YAML is a human-readable data serialization standard that is commonly used for configuration files. The major fields can be specified in numbers, floats, strings, lists and dictionaries.

In this example, the Pyharp module has one species, H2mix, and the atomic composition of it is specified as a dictionary with the key as the name of the atom and the value as its abundance. We use this entry to calculate its molecular weight.

Next, we specify the opacity sources, which are dictionary entries with the name of the opacity source as the key and its properties as the value. In this case, we have two opacity sources:

1. H2-molecule, which is a pyharp.opacity.MultiBand molecular opacity source with the data file sonora_2020_feh+000_co_100.data.196.pt and the dependent species is H2mix. We use these information to retrieve the opacity data from the file and calculate species indices in a tensor data.

2. H2-continuum, which is a pyharp.opacity.WaveTemp continuum opacity source with the data files H2-H2-eq.xiz.pt and H2-He-eq.xiz.pt and the dependent species is also H2mix. The data files for continuum absorption are shipped with Pyharp package. The following statement loads the data file directory such that Pyharp can find them:

from pyharp import h2_cia_legacy

Next, we specify the radiation bands, which is a list of band names. For each band, a dictionary entry is created with the band name as the key and its properties as the value.

Each band can have its own spectral range, opacity sources, solver, integration method and flags associated with radiative transfer computation. In this example, we have one band, sonora196, which has a spectral range of 30.8 to 38300 cm \(^{-1}\); the opacity sources are H2-molecule and H2-continuum; the radiative transfer solver is pydisort.Disort; the integration method is weight (correlated-K); and the flags passed to pydisort.Disort are lamber, quiet, onlyfl and planck, which assumes lambertian surface, quiet mode, flux-only calculation and activate planck function for the radiation source function respectively.

Please see the pydisort documentation for more details on the available options of pydisort.DisortOptions and their meanings.

Load modules

Then, we create a python file and load system modules include torch, numpy and os. These provide basic data structures and functions for numerical computation and file handling.

import torch
import os
import numpy as np
from typing import Tuple

We will download the pyharp.sonora opacity dataset and preprocess it to a format that can be used by the Pyharp module. These are the submodule functions we need.

from pyharp.sonora import (
        load_sonora_data,
        load_sonora_window,
        save_sonora_multiband,
        )

Last, we load Pyharp classes and functions. Pyharp uses a minimal set of classes and functions that are used to configure the radiation model and run the radiative transfer calculation.

import pyharp
from pyharp import (
        h2_cia_legacy,
        constants,
        RadiationOptions,
        Radiation,
        calc_dz_hypsometric,
        )

That’s all you need for the header of the program.

Pre-process sonora2020 data

The pyharp.sonora data is a large dataset that contains the opacity data for hydrogen and helium atmospheres. We first download a sample of it via the script

fetch-sonora --feh=+000 --co=100

The fetch-sonora script is a command line tool shipped with Pyharp that downloads the pyharp.sonora data from the Zenodo archive. This command will download a file named sonora_2020_feh+000_co_100.data.196.tar.gz in the current directory. The feh argument specifies the metallicity of the atmosphere and the co argument specifies the carbon-to-oxygen ratio (solar). You can check out the various options of the fetch-sonora script by running

fetch-sonora -h

def preprocess_sonora(fname: str):
    # returns dictionary of data
    data = load_sonora_data(fname)

    # save to ".pt" file
    save_sonora_multiband(fname, data, clean=False)

The above code block transforms the downloaded data into a .pt file, which is a binary format used by PyTorch for storing tensors and models.

Construct atmosphere

The following function constructs a simple atmosphere model with a given pressure and temperature point in the troposphere. It constructs an adiabatic atmosphere with a temperature profile that decreases with height, and an isothermal atmosphere with a constant temperature at a specified value.

def construct_atm(pmax: float, pmin: float,
                  ncol: int = 1,
                  nlyr: int = 100) -> dict[str, torch.Tensor]:
    p1bar, T1bar, Tmin = 1.e5, 169., 135.
    pres = torch.logspace(np.log10(pmax), np.log10(pmin), nlyr + 1, dtype=torch.float64)

    # adibatic
    temp = T1bar * torch.pow(pres / p1bar, 2. / 7.)

    # isothermal
    temp.clip_(min=Tmin)

    atm = {
        'pres' : (pres[1:] * pres[:-1]).sqrt().unsqueeze(0).expand(ncol, nlyr),
        'temp' : (temp[1:] * temp[:-1]).sqrt().unsqueeze(0).expand(ncol, nlyr),
    }
    return atm

Configure radiation bands

This function is where we configure the entire radiation model. Pyharp allows for multiple bands to be configured in a single radiation model. So, generally, we need to loop over the bands and configure them one by one.

There is no rule for dividing or defining the bands. However, each band must have the same number of spectral points, absorbers, radiation flags and opacity sources. The radiative transfer calculation is computed in parallel for every spectral point within each band, but sequentially over bands. So, there is a performance penalty when the number of bands is too large.

Since the pyharp.sonora dataset is homogeneous despite that it contains opacities for 196 wavenumber bands, we can group all of them in a single band, sonora196, with two opacity sources, H2-molecule, that is of type pyharp.opacity.MultiBand and H2-continuum, that is of type pyharp.opacity.WaveTemp. This improves the computational efficiency as the 196 radiative transfer calculations are performed in parallel.

def configure_bands(config_file: str,
                    ncol: int = 1,
                    nlyr: int = 100,
                    nstr: int = 4) -> Radiation:
    rad_op = RadiationOptions.from_yaml(config_file)
    wmin, wmax = load_sonora_window()

    for band in rad_op.bands():
        if band.name() == "sonora196":
            # Query wavenumber and weight from the opacity source
            op = band.opacities()["H2-molecule"]
            band.wavenumber(list(op.query_wavenumber()))
            band.weight(list(op.query_weight()))
            nwave = len(band.wavenumber())
            ng = int(nwave / len(wmin))

            band.disort().accur(1.0e-4)

            data = [wmin] * ng
            band.disort().wave_lower([x for col in zip(*data) for x in col])

            data = [wmax] * ng
            band.disort().wave_upper([x for col in zip(*data) for x in col])
        else:
            raise ValueError(f"Unknown band: {band.name()}")

    return Radiation(rad_op)

Run radiative transfer

The last step is to run the radiative transfer model. The crucial part of the code is to define the boundary conditions, bc, of the radiative transfer equations. Specifically, we need to define the albedo, emissivity, temperature of the surfaces bounding the atmosphere. Some boundary conditions are band-dependent, which is reflected in the keys of the bc dictionary. Others are band-independent, such as the temperature of the surface and the top of the atmosphere.

Here comes the forward method of pyharp.Radiation, which takes the concentration (in mol m \(^{-3}\)), layer thickness (in m), boundary conditions and atmospheric state (pressure in pa and temperature in K) as inputs and returns the radiative fluxes.

Throughout the Pyharp code, we will be using SI units unless otherwise specified.

def run_rt(rad: Radiation, conc: torch.Tensor, dz: torch.Tensor,
           atm: dict[str, torch.Tensor]) -> Tuple[torch.Tensor, torch.Tensor, torch.Tensor]:
    ncol = conc.shape[0]
    bc = {}
    for band in rad.options.bands():
        nwave = len(band.wavenumber())
        bc[band.name() + "/albedo"] = torch.ones((nwave, ncol), dtype=torch.float64)
        bc[band.name() + "/temis"] = torch.zeros((nwave, ncol), dtype=torch.float64)

    bc["btemp"] = torch.zeros((ncol), dtype=torch.float64)
    bc["ttemp"] = torch.zeros((ncol), dtype=torch.float64)

    return rad.forward(conc, dz, bc, atm)

The main program

This is the main program that calls the above functions and runs the radiative transfer model. The run_rt function produces three radiative fluxes:

1. net flux at each atmospheric layer

2. downward flux to the surface

3. upward flux at the top of the atmosphere

The first two fluxes are important for forcing the atmospheric circulation and the coupling between the atmosphere and the surface, while the last one is important for the energy balance of the planet.

if __name__ == "__main__":
    # prepare sonora2020 opacity data
    fname = "sonora_2020_feh+000_co_100.data.196"
    if not os.path.exists(fname + ".pt"):
        preprocess_sonora(fname)

    # configure atmosphere model
    atm = construct_atm(100.e5, 10., ncol=1, nlyr=100)

    # configure radiation model
    config_file = "example_sonora_2020.yaml"
    rad = configure_bands(config_file, ncol=1,
                          nlyr=atm['pres'].shape[-1], nstr=8)
    print(rad.options)

    # calculate concentration and layer thickness
    mean_mol_weight = pyharp.species_weights()[0]
    print("mean mol weight = ", mean_mol_weight)
    grav = 24.8 # m/s^2
    dz = calc_dz_hypsometric(atm["pres"], atm["temp"],
                             torch.tensor(mean_mol_weight * grav / constants.Rgas)
                             )
    print("dz = ", dz)
    conc = atm["pres"] / (atm["temp"] * constants.Rgas)
    conc.unsqueeze_(-1)

    # run rt
    netflux, dnflux, upflux = run_rt(rad, conc, dz, atm)
    print("netflux = ", netflux)
    print("surface flux = ", dnflux)
    print("toa flux = ", upflux)

The complete python code can be downloaded from this link and the yaml configuration file from this link.

References

[1]

Li, C., Le, T., Zhang, X., & Yung, Y. L. (2018). A high-performance atmospheric radiation package: With applications to the radiative energy budgets of giant planets. Journal of Quantitative Spectroscopy and Radiative Transfer, 217, 353-362.