vulcan.op

`op`

Module that includes all the numerical functions of VULCAN.

ReadRate() reads in the chemical network and construct the rate constants based on the T-P structure.
Integration() is the backbone of integrating for one time step
ODESolver() contains the functions for solving system of ODEs (e.g. dy/dt, Jacobian, etc.)

Imports

vulcan.config: Config class for handling configuration settings.
vulcan.paths: Paths to various files and directories used in VULCAN.
vulcan.build_atm: compo and compo_row functions for building the atmospheric composition.
vulcan.phy_const: Physical constants used in VULCAN.
vulcan.chem_funs: Functions related to the chemical network.

`ReadRate(vulcan_cfg)`

Bases: object

to read in rate constants from the network file and compute the reaction rates for the corresponding Tco and pco

Source code in src/vulcan/op.py

def __init__(self, vulcan_cfg: Config):

    self.cfg = vulcan_cfg
    self.i = 1
    # flag of trimolecular reaction
    self.re_tri, self.re_tri_k0 = False, False
    self.list_tri = []

`cfg = vulcan_cfg` `instance-attribute`

`i = 1` `instance-attribute`

`list_tri = []` `instance-attribute`

`lim_lowT_rates(var, atm)`

Source code in src/vulcan/op.py

def lim_lowT_rates(
    self, var, atm
):  # for setting up the lower limit of rate coefficients for low T
    for i in range(1, nr, 2):
        if var.Rf[i] == 'H + CH3 + M -> CH4 + M':
            T_mask = atm.Tco <= 277.5
            k0 = 6e-29
            kinf = 2.06e-10 * atm.Tco**-0.4  # from Moses+2005
            lowT_lim = k0 / (1.0 + k0 * atm.M / kinf)
            log.debug('using the low temperature limit for CH3 + H + M -> CH4 + M')
            log.debug('capping ')
            log.debug(var.k[i][T_mask])
            log.debug('at ')
            log.debug(lowT_lim[T_mask])
            var.k[i][T_mask] = lowT_lim[T_mask]

        elif var.Rf[i] == 'H + C2H4 + M -> C2H5 + M':
            T_mask = atm.Tco <= 300
            log.debug('using the low temperature limit for H + C2H4 + M -> C2H5 + M')
            log.debug('capping ')
            log.debug(var.k[i][T_mask])
            log.debug('at ')
            log.debug(3.7e-30)
            var.k[i][T_mask] = 3.7e-30  # from Moses+2005

        elif var.Rf[i] == 'H + C2H5 + M -> C2H6 + M':
            T_mask = atm.Tco <= 200
            log.debug('using the low temperature limit for H + C2H5 + M -> C2H6 + M')
            log.debug('capping ')
            log.debug(var.k[i][T_mask])
            log.debug('at ')
            log.debug(2.49e-27)
            var.k[i][T_mask] = 2.49e-27  # from Moses+2005

    return var

`make_bins_read_cross(var, atm)`

determining the bin range and only use the min and max wavelength that the molecules absorb to avoid photons with w0=1 (pure scatteing) in certain wavelengths var.cross stores the total absorption cross sections of each species, e.g. var.cross['H2O'] var.cross stores the IDIVIDUAL photodissociation cross sections for each bracnh, e.g. var.cross_J[('H2O',1)], which is equvilent to var.cross['H2O'] times the branching ratio of branch 1

Source code in src/vulcan/op.py

def make_bins_read_cross(self, var, atm):
    """
    determining the bin range and only use the min and max wavelength that the molecules absorb
    to avoid photons with w0=1 (pure scatteing) in certain wavelengths
    var.cross stores the total absorption cross sections of each species, e.g. var.cross['H2O']
    var.cross stores the IDIVIDUAL photodissociation cross sections for each bracnh, e.g. var.cross_J[('H2O',1)], which is equvilent to var.cross['H2O'] times the branching ratio of branch 1
    """
    photo_sp = list(var.photo_sp)
    ion_sp = list(var.ion_sp)
    absp_sp = photo_sp + ion_sp
    sp0 = photo_sp[0]
    nz = self.cfg.nz

    cross_raw, scat_raw = {}, {}
    ratio_raw, ion_ratio_raw = {}, {}
    cross_T_raw = {}

    # In the end, we do not need photons beyond the longest-wavelength threshold from all species (different from absorption)
    sp_label = np.genfromtxt(
        CROSS_DIR + 'thresholds.txt', dtype=str, usecols=0
    )  # taking the first column as labels
    lmd_data = np.genfromtxt(CROSS_DIR + 'thresholds.txt', skip_header=1)[
        :, 1
    ]  # discarding the fist column

    # for setting up the wavelength coverage
    threshold = {
        label: row for label, row in zip(sp_label, lmd_data) if label in species
    }  # only include the species in the current network
    var.threshold = threshold

    # reading in cross sections into dictionaries
    for n, sp in enumerate(absp_sp):
        if self.cfg.use_ion:
            try:
                cross_raw[sp] = np.genfromtxt(
                    CROSS_DIR + sp + '/' + sp + '_cross.csv',
                    dtype=float,
                    delimiter=',',
                    skip_header=1,
                    names=['lambda', 'cross', 'disso', 'ion'],
                )
            except FileNotFoundError:
                raise RuntimeError('Missing the cross section from ' + sp)
            if sp in ion_sp:
                try:
                    ion_ratio_raw[sp] = np.genfromtxt(
                        CROSS_DIR + sp + '/' + sp + '_ion_branch.csv',
                        dtype=float,
                        delimiter=',',
                        skip_header=1,
                        names=True,
                    )
                except FileNotFoundError:
                    raise RuntimeError('Missing the ion branching ratio from ' + sp)
        else:
            try:
                cross_raw[sp] = np.genfromtxt(
                    CROSS_DIR + sp + '/' + sp + '_cross.csv',
                    dtype=float,
                    delimiter=',',
                    skip_header=1,
                    names=['lambda', 'cross', 'disso'],
                )
            except FileNotFoundError:
                raise RuntimeError('Missing the cross section from ' + sp)

        # reading in the branching ratios
        # for i in range(1,var.n_branch[sp]+1): # branch index should start from 1
        if sp in photo_sp:  # excluding ion_sp
            try:
                ratio_raw[sp] = np.genfromtxt(
                    CROSS_DIR + sp + '/' + sp + '_branch.csv',
                    dtype=float,
                    delimiter=',',
                    skip_header=1,
                    names=True,
                )
            except FileNotFoundError:
                raise RuntimeError('Missing the branching ratio from ' + sp)

        # reading in temperature dependent cross sections
        if sp in self.cfg.T_cross_sp:
            T_list = []
            for temp_file in os.listdir('thermo/photo_cross/' + sp + '/'):
                if temp_file.startswith(sp) and temp_file.endswith('K.csv'):
                    temp = temp_file
                    temp = temp.replace(sp, '')
                    temp = temp.replace('_cross_', '')
                    temp = temp.replace('K.csv', '')
                    T_list.append(int(temp))
                    var.cross_T_sp_list[sp] = T_list
            for tt in T_list:
                if self.cfg.use_ion:  # usually the T-dependent cross sections are only measured in the photodissociation-relavent wavelengths so cross_tot = cross_diss
                    cross_T_raw[(sp, tt)] = np.genfromtxt(
                        CROSS_DIR + sp + '/' + sp + '_cross_' + str(tt) + 'K.csv',
                        dtype=float,
                        delimiter=',',
                        skip_header=1,
                        names=['lambda', 'cross', 'disso', 'ion'],
                    )
                else:
                    cross_T_raw[(sp, tt)] = np.genfromtxt(
                        CROSS_DIR + sp + '/' + sp + '_cross_' + str(tt) + 'K.csv',
                        dtype=float,
                        delimiter=',',
                        skip_header=1,
                        names=['lambda', 'cross', 'disso'],
                    )
            # room-T cross section
            cross_T_raw[(sp, 300)] = cross_raw[sp]
            var.cross_T_sp_list[sp].append(300)

        if cross_raw[sp]['cross'][0] == 0 or cross_raw[sp]['cross'][-1] == 0:
            raise IOError('Please remove the zeros in the cross file of ' + sp)

        if n == 0:  # the first species
            bin_min = cross_raw[sp]['lambda'][0]
            bin_max = cross_raw[sp]['lambda'][-1]
            # photolysis threshold
            try:
                diss_max = threshold[sp]
            except (KeyError, IndexError):
                raise RuntimeError(sp + ' not in threshol.txt')

        else:
            sp_min, sp_max = cross_raw[sp]['lambda'][0], cross_raw[sp]['lambda'][-1]
            if sp_min < bin_min:
                bin_min = sp_min
            if sp_max > bin_max:
                bin_max = sp_max
            try:
                if threshold[sp] > diss_max:
                    diss_max = threshold[sp]
            except (KeyError, IndexError):
                raise RuntimeError(sp + ' not in threshol.txt')

    # constraining the bin_min and bin_max by the default values defined in store.py
    bin_min = max(bin_min, var.def_bin_min)
    bin_max = min(bin_max, var.def_bin_max, diss_max)
    log.info(
        'Input stellar spectrum from '
        + '{:.1f}'.format(var.def_bin_min)
        + ' to '
        + '{:.1f}'.format(var.def_bin_max)
    )
    log.debug('Photodissociation threshold: ' + '{:.1f}'.format(diss_max))
    log.info(
        'Using wavelength bins from ' + '{:.1f}'.format(bin_min) + ' to ' + str(bin_max)
    )

    var.dbin1 = self.cfg.dbin1
    var.dbin2 = self.cfg.dbin2
    if self.cfg.dbin_12trans >= bin_min and self.cfg.dbin_12trans <= bin_max:
        bins = np.concatenate(
            (
                np.arange(bin_min, self.cfg.dbin_12trans, var.dbin1),
                np.arange(self.cfg.dbin_12trans, bin_max, var.dbin2),
            )
        )
    else:
        bins = np.arange(bin_min, bin_max, var.dbin1)
    var.bins = bins
    var.nbin = len(bins)

    # all variables that depend on the size of nbins
    # the direct beam (staggered)
    var.sflux = np.zeros((nz + 1, var.nbin))
    # the diffusive flux (staggered)
    var.dflux_u, var.dflux_d = np.zeros((nz + 1, var.nbin)), np.zeros((nz + 1, var.nbin))
    # the total actinic flux (non-staggered)
    var.aflux = np.zeros((nz, var.nbin))
    # the total actinic flux from the previous calculation
    prev_aflux = np.zeros((nz, var.nbin))

    # staggered
    var.tau = np.zeros((nz + 1, var.nbin))
    # the stellar flux at TOA
    var.sflux_top = np.zeros(var.nbin)

    # read_cross
    # creat a dict of cross section with key=sp and values=bins in data_var
    var.cross = dict(
        [(sp, np.zeros(var.nbin)) for sp in absp_sp]
    )  # including photo_sp and ion_sp

    # read cross of disscoiation
    var.cross_J = dict(
        [
            ((sp, i), np.zeros(var.nbin))
            for sp in photo_sp
            for i in range(1, var.n_branch[sp] + 1)
        ]
    )
    var.cross_scat = dict([(sp, np.zeros(var.nbin)) for sp in self.cfg.scat_sp])

    # for temperature-dependent cross sections
    var.cross_T = dict([(sp, np.zeros((nz, var.nbin))) for sp in self.cfg.T_cross_sp])
    var.cross_J_T = dict(
        [
            ((sp, i), np.zeros((nz, var.nbin)))
            for sp in self.cfg.T_cross_sp
            for i in range(1, var.n_branch[sp] + 1)
        ]
    )

    # read cross of ionisation
    if self.cfg.use_ion:
        var.cross_Jion = dict(
            [
                ((sp, i), np.zeros(var.nbin))
                for sp in ion_sp
                for i in range(1, var.ion_branch[sp] + 1)
            ]
        )

    for sp in (
        photo_sp
    ):  # photodissociation only; photoionization takes a separate branch ratio file
        # for values outside the boundary => fill_value = 0
        inter_cross = interpolate.interp1d(
            cross_raw[sp]['lambda'],
            cross_raw[sp]['cross'],
            bounds_error=False,
            fill_value=0,
        )
        inter_cross_J = interpolate.interp1d(
            cross_raw[sp]['lambda'],
            cross_raw[sp]['disso'],
            bounds_error=False,
            fill_value=0,
        )
        inter_ratio = {}  # excluding ionization branches

        for i in range(
            1, var.n_branch[sp] + 1
        ):  # fill_value extends the first and last elements for branching ratios
            br_key = 'br_ratio_' + str(i)
            try:
                inter_ratio[i] = interpolate.interp1d(
                    ratio_raw[sp]['lambda'],
                    ratio_raw[sp][br_key],
                    bounds_error=False,
                    fill_value=(ratio_raw[sp][br_key][0], ratio_raw[sp][br_key][-1]),
                )
            except (ValueError, KeyError, IndexError):
                log.error(
                    'The branches in the network file does not match the branchong ratio file for '
                    + str(sp)
                )

        # using a loop instead of an array because it's easier to handle the branching ratios
        for n, ld in enumerate(bins):
            var.cross[sp][n] = inter_cross(ld)

            # using the branching ratio (from the files) to construct the individual cross section of each branch
            for i in range(1, var.n_branch[sp] + 1):
                var.cross_J[(sp, i)][n] = inter_cross_J(ld) * inter_ratio[i](ld)

        # make var.cross_T[(sp,i)] and var.cross_J_T[(sp,i)] here in 2D array: nz * bins (same shape as tau)
        # T-dependent cross sections are usually only measured in the photodissociation-relavent wavelengths so cross_tot = cross_diss
        if sp in self.cfg.T_cross_sp:
            # T list of species sp that have T-depedent cross sections (inclduing 300 K for inter_cross)
            T_list = np.array(var.cross_T_sp_list[sp])
            max_T_sp = np.amax(T_list)
            min_T_sp = np.amin(T_list)

            for lev, Tz in enumerate(atm.Tco):  # looping z
                Tz_between = False  # flag for Tz in between any two elements in T_list
                # define the interpolating T range
                if list(T_list[T_list <= Tz]) and list(T_list[T_list > Tz]):
                    Tlow = T_list[T_list <= Tz].max()  # closest T in T_list smaller than Tz
                    Thigh = T_list[T_list > Tz].min()  # closest T in T_list larger than Tz
                    Tz_between = True

                    # find the wavelength range that are included in both cross_T_raw[(sp,Tlow)] and cross_T_raw[(sp,Thigh)]
                    ld_min = max(
                        cross_T_raw[(sp, Tlow)]['lambda'][0],
                        cross_T_raw[(sp, Thigh)]['lambda'][0],
                    )
                    ld_max = min(
                        cross_T_raw[(sp, Tlow)]['lambda'][-1],
                        cross_T_raw[(sp, Thigh)]['lambda'][-1],
                    )
                    inter_cross_lowT = interpolate.interp1d(
                        cross_T_raw[(sp, Tlow)]['lambda'],
                        cross_T_raw[(sp, Tlow)]['cross'],
                        bounds_error=False,
                        fill_value=0,
                    )
                    inter_cross_highT = interpolate.interp1d(
                        cross_T_raw[(sp, Thigh)]['lambda'],
                        cross_T_raw[(sp, Thigh)]['cross'],
                        bounds_error=False,
                        fill_value=0,
                    )
                    inter_cross_J_lowT = interpolate.interp1d(
                        cross_T_raw[(sp, Tlow)]['lambda'],
                        cross_T_raw[(sp, Tlow)]['disso'],
                        bounds_error=False,
                        fill_value=0,
                    )
                    inter_cross_J_highT = interpolate.interp1d(
                        cross_T_raw[(sp, Thigh)]['lambda'],
                        cross_T_raw[(sp, Thigh)]['disso'],
                        bounds_error=False,
                        fill_value=0,
                    )

                    for n, ld in enumerate(bins):  # looping bins
                        # not within the T-cross wavelength range
                        if ld < ld_min or ld > ld_max:
                            var.cross_T[sp][lev, n] = var.cross[sp][n]
                            # don't forget the cross_J_T branches
                            for i in range(1, var.n_branch[sp] + 1):
                                var.cross_J_T[(sp, i)][lev, n] = var.cross_J[(sp, i)][n]

                        else:
                            # update: inerpolation in log10 for cross sections and linearly between Tlow and Thigh
                            log_lowT = np.log10(inter_cross_lowT(ld))
                            log_highT = np.log10(inter_cross_highT(ld))
                            if np.isinf(log_lowT):
                                log_lowT = -100.0  # replacing -inf with -100
                            if np.isinf(log_highT):
                                log_highT = -100.0

                            inter_T = interpolate.interp1d(
                                [Tlow, Thigh], [log_lowT, log_highT], axis=0
                            )  # at wavelength ld, interpolating between Tlow and Thigh in log10
                            if inter_T(Tz) == -100:
                                var.cross_T[sp][lev, n] == 0.0
                            else:
                                var.cross_T[sp][lev, n] = 10 ** (inter_T(Tz))

                            # update: inerpolation in log10 for cross sections and linearly between Tlow and Thigh
                            # using the branching ratio (from the files) to construct the individual cross section of each branch
                            for i in range(1, var.n_branch[sp] + 1):
                                J_log_lowT = np.log10(inter_cross_J_lowT(ld))
                                J_log_highT = np.log10(inter_cross_J_highT(ld))
                                if np.isinf(J_log_lowT):
                                    J_log_lowT = -100.0  # replacing -inf with -100
                                if np.isinf(J_log_highT):
                                    J_log_highT = -100.0

                                inter_cross_J_T = interpolate.interp1d(
                                    [Tlow, Thigh], [J_log_lowT, J_log_highT], axis=0
                                )

                                if inter_cross_J_T(Tz) == -100:
                                    var.cross_J_T[(sp, i)][lev, n] = 0.0
                                else:
                                    var.cross_J_T[(sp, i)][lev, n] = 10 ** (
                                        inter_cross_J_T(Tz)
                                    ) * inter_ratio[i](
                                        ld
                                    )  # same inter_ratio[i](ld) as the standard one above

                elif not list(
                    T_list[T_list < Tz]
                ):  # Tz equal or smaller than all T in T_list including 300K (empty list)
                    if min_T_sp == 300:
                        var.cross_T[sp][lev] = var.cross[
                            sp
                        ]  # using the cross section at room T
                        for i in range(1, var.n_branch[sp] + 1):
                            var.cross_J_T[(sp, i)][lev] = var.cross_J[(sp, i)]
                    else:  # min_T_sp != 300; T-cross lower than room temperature
                        # the wavelength range of cross_T_raw at T = min_T_sp
                        ld_min, ld_max = (
                            cross_T_raw[(sp, min_T_sp)]['lambda'][0],
                            cross_T_raw[(sp, min_T_sp)]['lambda'][-1],
                        )
                        inter_cross_lowT = interpolate.interp1d(
                            cross_T_raw[(sp, min_T_sp)]['lambda'],
                            cross_T_raw[(sp, min_T_sp)]['cross'],
                            bounds_error=False,
                            fill_value=0,
                        )
                        inter_cross_J_lowT = interpolate.interp1d(
                            cross_T_raw[(sp, min_T_sp)]['lambda'],
                            cross_T_raw[(sp, min_T_sp)]['disso'],
                            bounds_error=False,
                            fill_value=0,
                        )
                        for n, ld in enumerate(bins):  # looping bins
                            # not within the T-cross wavelength range
                            if ld < ld_min or ld > ld_max:
                                var.cross_T[sp][lev, n] = var.cross[sp][n]
                                # don't forget the cross_J_T branches
                                for i in range(1, var.n_branch[sp] + 1):
                                    var.cross_J_T[(sp, i)][lev, n] = var.cross_J[(sp, i)][n]
                            else:
                                var.cross_T[sp][lev, n] = inter_cross_lowT(ld)
                                # using the branching ratio (from the files) to construct the individual cross section of each branch
                                for i in range(1, var.n_branch[sp] + 1):
                                    var.cross_J_T[(sp, i)][lev, n] = inter_cross_J_lowT(
                                        ld
                                    ) * inter_ratio[i](
                                        ld
                                    )  # same inter_ratio[i](ld) as the standard one above

                else:  # Tz equal or larger than all T in T_list (empty list)
                    # the wavelength range of cross_T_raw[(sp,Thigh)]

                    if max_T_sp == 300:
                        var.cross_T[sp][lev] = var.cross[
                            sp
                        ]  # using the cross section at room T
                        for i in range(1, var.n_branch[sp] + 1):
                            var.cross_J_T[(sp, i)][lev] = var.cross_J[(sp, i)]
                    else:  # the wavelength range of cross_T_raw at T = max_T_sp
                        ld_min, ld_max = (
                            cross_T_raw[(sp, max_T_sp)]['lambda'][0],
                            cross_T_raw[(sp, max_T_sp)]['lambda'][-1],
                        )
                        inter_cross_highT = interpolate.interp1d(
                            cross_T_raw[(sp, max_T_sp)]['lambda'],
                            cross_T_raw[(sp, max_T_sp)]['cross'],
                            bounds_error=False,
                            fill_value=0,
                        )
                        inter_cross_J_highT = interpolate.interp1d(
                            cross_T_raw[(sp, max_T_sp)]['lambda'],
                            cross_T_raw[(sp, max_T_sp)]['disso'],
                            bounds_error=False,
                            fill_value=0,
                        )
                        for n, ld in enumerate(bins):  # looping bins
                            # not within the T-cross wavelength range
                            if ld < ld_min or ld > ld_max:
                                var.cross_T[sp][lev, n] = var.cross[sp][n]
                                # don't forget the cross_J_T branches
                                for i in range(1, var.n_branch[sp] + 1):
                                    var.cross_J_T[(sp, i)][lev, n] = var.cross_J[(sp, i)][n]
                            else:
                                var.cross_T[sp][lev, n] = inter_cross_highT(ld)

                                # using the branching ratio (from the files) to construct the individual cross section of each branch
                                for i in range(1, var.n_branch[sp] + 1):
                                    var.cross_J_T[(sp, i)][lev, n] = inter_cross_J_highT(
                                        ld
                                    ) * inter_ratio[i](
                                        ld
                                    )  # same inter_ratio[i](ld) as the standard one above

    if self.cfg.use_ion:
        for sp in ion_sp:
            if sp not in photo_sp:
                inter_cross = interpolate.interp1d(
                    cross_raw[sp]['lambda'],
                    cross_raw[sp]['cross'],
                    bounds_error=False,
                    fill_value=0,
                )

            inter_cross_Jion = interpolate.interp1d(
                cross_raw[sp]['lambda'],
                cross_raw[sp]['ion'],
                bounds_error=False,
                fill_value=0,
            )
            ion_inter_ratio = {}  # For ionization branches

            for i in range(
                1, var.ion_branch[sp] + 1
            ):  # fill_value extends the first and last elements for branching ratios
                br_key = 'br_ratio_' + str(i)
                try:
                    ion_inter_ratio[i] = interpolate.interp1d(
                        ion_ratio_raw[sp]['lambda'],
                        ion_ratio_raw[sp][br_key],
                        bounds_error=False,
                        fill_value=(
                            ion_ratio_raw[sp][br_key][0],
                            ion_ratio_raw[sp][br_key][-1],
                        ),
                    )
                except (KeyError, IndexError, ValueError):
                    log.error(
                        'The ionic branches in the network file does not match the branchong ratio file for '
                        + str(sp)
                    )

            for n, ld in enumerate(bins):
                # for species noe appeared in photodissociation but only in photoionization, like H
                if sp not in photo_sp:
                    var.cross[sp][n] = inter_cross(ld)
                for i in range(1, var.ion_branch[sp] + 1):
                    var.cross_Jion[(sp, i)][n] = inter_cross_Jion(ld) * ion_inter_ratio[i](
                        ld
                    )
    # end of if self.cfg.use_ion :

    # reading in cross sections of Rayleigh Scattering
    for sp in self.cfg.scat_sp:
        scat_raw[sp] = np.genfromtxt(
            CROSS_DIR + 'rayleigh/' + sp + '_scat.txt',
            dtype=float,
            skip_header=1,
            names=['lambda', 'cross'],
        )

        # for values outside the boundary => fill_value = 0
        inter_scat = interpolate.interp1d(
            scat_raw[sp]['lambda'], scat_raw[sp]['cross'], bounds_error=False, fill_value=0
        )

        for n, ld in enumerate(bins):
            var.cross_scat[sp][n] = inter_scat(ld)

`read_rate(var, atm)`

Source code in src/vulcan/op.py

def read_rate(self, var, atm):

    nz = self.cfg.nz

    (
        Rf,
        Rindx,
        a,
        n,
        E,
        a_inf,
        n_inf,
        E_inf,
        k,
        k_fun,
        k_inf,
        kinf_fun,
        k_fun_new,
        pho_rate_index,
    ) = (
        var.Rf,
        var.Rindx,
        var.a,
        var.n,
        var.E,
        var.a_inf,
        var.n_inf,
        var.E_inf,
        var.k,
        var.k_fun,
        var.k_inf,
        var.kinf_fun,
        var.k_fun_new,
        var.pho_rate_index,
    )
    ion_rate_index = var.ion_rate_index

    i = self.i
    re_tri, re_tri_k0 = self.re_tri, self.re_tri_k0
    list_tri = self.list_tri

    Tco = atm.Tco
    M = atm.M.copy()

    special_re = False
    conden_re = False
    recomb_re = False
    photo_re = False
    ion_re = False
    # end_re = False
    # br_read  = False
    # ion_br_read = False

    photo_sp = []
    ion_sp = []

    with open(self.cfg.network) as f:
        all_lines = f.readlines()
        for line_indx, line in enumerate(all_lines):
            # switch to 3-body and dissociation reations
            if line.startswith('# 3-body'):
                re_tri = True

            if line.startswith('# 3-body reactions without high-pressure rates'):
                re_tri_k0 = True

            elif line.startswith('# special'):
                re_tri = False
                re_tri_k0 = False
                special_re = True  # switch to reactions with special forms (hard coded)

            elif line.startswith('# condensation'):
                re_tri = False
                re_tri_k0 = False
                special_re = False
                conden_re = True
                var.conden_indx = i

            elif line.startswith('# radiative'):
                re_tri = False
                re_tri_k0 = False
                special_re = False
                conden_re = False
                recomb_re = True
                var.recomb_indx = i

            elif line.startswith('# photo'):
                re_tri = False
                re_tri_k0 = False
                special_re = False  # turn off reading in the special form
                conden_re = False
                recomb_re = False
                photo_re = True
                var.photo_indx = i

            elif line.startswith('# ionisation'):
                re_tri = False
                re_tri_k0 = False
                special_re = False  # turn off reading in the special form
                conden_re = False
                recomb_re = False
                photo_re = False
                ion_re = True
                var.ion_indx = i

            elif line.startswith('# reverse stops'):
                var.special_re = False
                var.stop_rev_indx = i

            # skip common lines and blank lines
            # ========================================================================================
            if (
                not line.startswith('#')
                and line.strip()
                and not special_re
                and not conden_re
                and not photo_re
                and not ion_re
            ):  # if not starts
                Rf[i] = line.partition('[')[-1].rpartition(']')[0].strip()
                li = line.partition(']')[-1].strip()
                columns = li.split()
                Rindx[i] = int(line.partition('[')[0].strip())
                a[i] = float(columns[0])
                n[i] = float(columns[1])
                E[i] = float(columns[2])

                # switching to trimolecular reactions (len(columns) > 3 for those with high-P limit rates)
                if re_tri and not re_tri_k0:
                    a_inf[i] = float(columns[3])
                    n_inf[i] = float(columns[4])
                    E_inf[i] = float(columns[5])
                    list_tri.append(i)

                if columns[-1].strip() == 'He':
                    re_He = i
                elif columns[-1].strip() == 'ex1':
                    re_CH3OH = i

                # Note: make the defaut i=i
                k_fun[i] = lambda temp, mm, i=i: a[i] * temp ** n[i] * np.exp(-E[i] / temp)

                if not re_tri:
                    k[i] = k_fun[i](Tco, M)

                # for 3-body reactions, also calculating k_inf
                elif re_tri and len(columns) >= 6:
                    kinf_fun[i] = lambda temp, i=i: (
                        a_inf[i] * temp ** n_inf[i] * np.exp(-E_inf[i] / temp)
                    )
                    k_fun_new[i] = lambda temp, mm, i=i: (
                        (a[i] * temp ** n[i] * np.exp(-E[i] / temp))
                        / (
                            1
                            + (a[i] * temp ** n[i] * np.exp(-E[i] / temp))
                            * mm
                            / (a_inf[i] * temp ** n_inf[i] * np.exp(-E_inf[i] / temp))
                        )
                    )

                    # k[i] = k_fun_new[i](Tco, M)
                    k_inf = a_inf[i] * Tco ** n_inf[i] * np.exp(-E_inf[i] / Tco)
                    k[i] = k_fun[i](Tco, M)
                    k[i] = k[i] / (1 + k[i] * M / k_inf)

                else:  # for 3-body reactions without high-pressure rates
                    k[i] = k_fun[i](Tco, M)

                i += 2
                # end if not
            # ========================================================================================
            elif special_re and line.strip() and not line.startswith('#'):
                Rindx[i] = int(line.partition('[')[0].strip())
                Rf[i] = line.partition('[')[-1].rpartition(']')[0].strip()

                if Rf[i] == 'OH + CH3 + M -> CH3OH + M':
                    log.debug('Using special form for the reaction: ' + Rf[i])

                    k[i] = 1.932e3 * Tco**-9.88 * np.exp(
                        -7544.0 / Tco
                    ) + 5.109e-11 * Tco**-6.25 * np.exp(-1433.0 / Tco)
                    k_inf = 1.031e-10 * Tco**-0.018 * np.exp(16.74 / Tco)
                    # the pressure dependence from Jasper 2017
                    Fc = (
                        0.1855 * np.exp(-Tco / 155.8)
                        + 0.8145 * np.exp(-Tco / 1675.0)
                        + np.exp(-4531.0 / Tco)
                    )
                    nn = 0.75 - 1.27 * np.log(Fc)
                    ff = np.exp(
                        np.log(Fc) / (1.0 + (np.log(k[i] * M / k_inf) / nn**2) ** 2)
                    )

                    k[i] = k[i] / (1 + k[i] * M / k_inf) * ff

                    k_fun[i] = lambda temp, mm, i=i: (
                        1.932e3 * temp**-9.88 * np.exp(-7544.0 / temp)
                        + 5.109e-11 * temp**-6.25 * np.exp(-1433.0 / temp)
                    )
                    kinf_fun[i] = lambda temp, mm, i=i: (
                        1.031e-10 * temp**-0.018 * np.exp(16.74 / temp)
                    )
                    k_fun_new[i] = lambda temp, mm, i=i: (
                        (
                            1.932e3 * temp**-9.88 * np.exp(-7544.0 / temp)
                            + 5.109e-11 * temp**-6.25 * np.exp(-1433.0 / temp)
                        )
                        / (
                            1
                            + (
                                1.932e3 * temp**-9.88 * np.exp(-7544.0 / temp)
                                + 5.109e-11 * temp**-6.25 * np.exp(-1433.0 / temp)
                            )
                            * mm
                            / (1.031e-10 * temp**-0.018 * np.exp(16.74 / temp))
                        )
                    )

                i += 2

            # Testing condensation
            elif conden_re and line.strip() and not line.startswith('#'):
                Rindx[i] = int(line.partition('[')[0].strip())
                Rf[i] = line.partition('[')[-1].rpartition(']')[0].strip()

                var.conden_re_list.append(i)
                k[i] = np.zeros(nz)
                k[i + 1] = np.zeros(nz)

                i += 2

            # setting photo dissociation reactions to zeros
            elif photo_re and line.strip() and not line.startswith('#'):
                k[i] = np.zeros(nz)
                Rf[i] = line.partition('[')[-1].rpartition(']')[0].strip()

                # adding the photo species
                photo_sp.append(Rf[i].split()[0])

                li = line.partition(']')[-1].strip()
                columns = li.split()
                Rindx[i] = int(line.partition('[')[0].strip())
                # columns[0]: the species being dissocited; branch index: columns[1]
                pho_rate_index[(columns[0], int(columns[1]))] = Rindx[i]

                # store the number of branches
                var.n_branch[columns[0]] = int(columns[1])

                i += 2

            # setting photo ionization reactions to zeros
            elif (
                ion_re and line.strip() and not line.startswith('#')
            ):  # and end_re == False
                k[i] = np.zeros(nz)
                Rf[i] = line.partition('[')[-1].rpartition(']')[0].strip()

                # chekcing if it already existed in the photo species
                ion_sp.append(Rf[i].split()[0])

                li = line.partition(']')[-1].strip()
                columns = li.split()
                Rindx[i] = int(line.partition('[')[0].strip())
                # columns[0]: the species being dissocited; branch index: columns[1]
                ion_rate_index[(columns[0], int(columns[1]))] = Rindx[i]

                # store the number of branches
                var.ion_branch[columns[0]] = int(columns[1])

                i += 2

    k_fun.update(k_fun_new)

    # store k into data_var
    # remeber k_fun has not removed reactions from remove_list
    var.k = k
    var.k_fun = k_fun
    var.kinf_fun = kinf_fun

    var.photo_sp = set(photo_sp)
    if self.cfg.use_ion:
        var.ion_sp = set(ion_sp)

    return var

`remove_rate(var)`

Source code in src/vulcan/op.py

def remove_rate(self, var):

    nz = self.cfg.nz

    for i in self.cfg.remove_list:
        var.k[i] = np.repeat(0.0, nz)
        var.k_fun[i] = lambda temp, mm, i=i: np.repeat(0.0, nz)

    return var

`rev_rate(var, atm)`

Source code in src/vulcan/op.py

def rev_rate(self, var, atm):

    nz = self.cfg.nz

    rev_list = range(2, var.stop_rev_indx, 2)
    # setting the rest reversal zeros
    for i in range(var.stop_rev_indx + 1, nr + 1, 2):
        var.k[i] = np.zeros(nz)

    Tco = atm.Tco

    # reversing rates and storing into data_var
    log.debug('Reverse rates from R1 to R' + str(var.stop_rev_indx - 2))
    log.debug('Rates greater than 1e-6:')
    for i in rev_list:
        if i in self.cfg.remove_list:
            var.k[i] = np.repeat(0.0, nz)
        else:
            var.k_fun[i] = lambda temp, mm, i=i: (
                var.k_fun[i - 1](temp, mm) / Gibbs(i - 1, temp)
            )
            var.k[i] = var.k[i - 1] / Gibbs(i - 1, Tco)

        if np.any(var.k[i] > 1.0e-6):
            log.debug('R' + str(i) + ' ' + var.Rf[i - 1] + ' :  ' + str(np.amax(var.k[i])))
        if np.any(var.k[i - 1] > 1.0e-6):
            log.debug(
                'R' + str(i - 1) + ' ' + var.Rf[i - 1] + ' :  ' + str(np.amax(var.k[i - 1]))
            )

    return var

`Integration(odesolver, output, vulcan_cfg)`

Bases: object

time-stepping until the stopping criteria (steady-state) is satisfied

all the operators required in the continuity equation: dn/dt + dphi/dz = P - L

or class incorporating the esential numerical operations?

Source code in src/vulcan/op.py

def __init__(self, odesolver, output, vulcan_cfg: Config):

    self.cfg = vulcan_cfg
    self.mtol = self.cfg.mtol
    self.atol = self.cfg.atol
    self.output = output

    self.odesolver = odesolver
    self.non_gas_sp = self.cfg.non_gas_sp
    self.use_settling = self.cfg.use_settling

    # warn about in-development features
    if self.cfg.use_vm_mol:
        log.warning(
            'New upwind scheme for molecular diffusion has been enabled in VULCAN (in development)'
        )
        raise RuntimeError(
            'You must disable the upwind scheme by setting `use_vm_mol=False`'
        )

    # import AGNI?
    if self.cfg.agni_call_frq > 0:
        from .agni import run_agni

        self.run_agni = run_agni

    # including photoionisation
    if self.cfg.use_photo:
        self.update_photo_frq = self.cfg.ini_update_photo_frq

    if self.cfg.use_condense:
        self.non_gas_sp_index = [species.index(sp) for sp in self.non_gas_sp]
        self.condense_sp_index = [species.index(sp) for sp in self.cfg.condense_sp]

`atol = self.cfg.atol` `instance-attribute`

`cfg = vulcan_cfg` `instance-attribute`

`condense_sp_index = [(species.index(sp)) for sp in (self.cfg.condense_sp)]` `instance-attribute`

`mtol = self.cfg.mtol` `instance-attribute`

`non_gas_sp = self.cfg.non_gas_sp` `instance-attribute`

`non_gas_sp_index = [(species.index(sp)) for sp in (self.non_gas_sp)]` `instance-attribute`

`odesolver = odesolver` `instance-attribute`

`output = output` `instance-attribute`

`run_agni = run_agni` `instance-attribute`

`update_photo_frq = self.cfg.ini_update_photo_frq` `instance-attribute`

`use_settling = self.cfg.use_settling` `instance-attribute`

`backup(var)`

Source code in src/vulcan/op.py

def backup(self, var):
    var.y_prev = np.copy(var.y)
    var.dy_prev = np.copy(var.dy)
    var.atom_loss_prev = var.atom_loss.copy()
    return var

`conden(var, atm)`

Updating the condensation reactions according to the new number density using the condensation growth timescale in the contiuum regime (not in the kinetic regime)

Note that when n_g -> n_s, n_s is still the number density of "molecules", not "particles." So I scaled down the evaporation rate by n_mol. n_s / n_mol should also be used for plotting.

Source code in src/vulcan/op.py

def conden(self, var, atm):
    """
    Updating the condensation reactions according to the new number density
    using the condensation growth timescale in the contiuum regime (not in the kinetic regime)

    Note that when n_g -> n_s, n_s is still the number density of "molecules", not "particles."
    So I scaled down the evaporation rate by n_mol.
    n_s / n_mol should also be used for plotting.
    """

    nz = self.cfg.nz

    for re in var.conden_re_list:
        if var.Rf[re] == 'H2O -> H2O_l_s' and 'H2O' in self.cfg.condense_sp:
            # using realxation for water condensation
            if self.cfg.use_relax:
                var.k[re] = np.repeat(0.0, nz)
                var.k[re + 1] = np.repeat(0.0, nz)
            else:
                m = 18.0 / Navo
                rho_p = atm.rho_p['H2O_l_s']
                r_p = atm.r_p['H2O_l_s']
                # relative humidity
                sat_humidity = atm.sat_p['H2O'] / kb / atm.Tco * self.cfg.humidity

                # this is based on the kinetic regime
                rate_c = (
                    m
                    / (4 * rho_p)
                    * (8 * kb * atm.Tco / np.pi / m) ** 0.5
                    * (var.y[:, species.index('H2O')] - sat_humidity)
                    / r_p
                )

                # new approach: contiuum regime DM/rho c
                Dg = np.insert(
                    atm.Dzz[:, species.index('H2O')], 0, atm.Dzz[0, species.index('H2O')]
                )
                rate = (
                    Dg
                    * m
                    / rho_p
                    / r_p**2
                    * (var.y[:, species.index('H2O')] - sat_humidity)
                )

                # how many gas molecules are contained in one particle with the assumed size r_p
                n_mol = 4.0 / 3 * np.pi * r_p**3 * rho_p / m

                var.k[re] = rate
                var.k[re + 1] = rate  # /n_mol

                # positive: condensation
                var.k[re] = np.maximum(var.k[re], 0)
                # negative: evaporation
                var.k[re + 1] = np.minimum(var.k[re + 1], 0)
                var.k[re + 1] = np.abs(var.k[re + 1])

        elif var.Rf[re] == 'NH3 -> NH3_l' and 'NH3' in self.cfg.condense_sp:
            # using realxation for water condensation
            if self.cfg.use_relax:
                var.k[re] = np.repeat(0.0, nz)
                var.k[re + 1] = np.repeat(0.0, nz)
            else:
                m = 17.0 / Navo
                rho_p = atm.rho_p['NH3_l_s']
                r_p = atm.r_p['NH3_l_s']  # assuming 1 micron

                # rate_c = m/(4*rho_p)*(8*kb*atm.Tco/np.pi/m)**0.5 *(var.y[:,species.index('NH3')]-atm.sat_p['NH3']/kb/atm.Tco)/r_p

                # new approach: contiuum regime DM/rho c
                Dg = np.insert(
                    atm.Dzz[:, species.index('NH3')], 0, atm.Dzz[0, species.index('NH3')]
                )
                rate = (
                    Dg
                    * m
                    / rho_p
                    / r_p**2
                    * (var.y[:, species.index('NH3')] - atm.sat_p['NH3'] / kb / atm.Tco)
                )

                # how many gas molecules are contained in one particle with the assumed size r_p
                n_mol = 4.0 / 3 * np.pi * r_p**3 * rho_p / m

                var.k[re] = rate
                var.k[re + 1] = rate  # /n_mol

                # positive: condensation
                var.k[re] = np.maximum(var.k[re], 0)
                # negative: evaporation
                var.k[re + 1] = np.minimum(var.k[re + 1], 0)
                var.k[re + 1] = np.abs(var.k[re + 1])

        elif var.Rf[re] == 'H2SO4 -> H2SO4_l' and 'H2SO4' in self.cfg.condense_sp:
            m = 98.022 / Navo
            rho_p = atm.rho_p['H2SO4_l']
            r_p = atm.r_p['H2SO4_l']

            # new approach: contiuum regime DM/rho c
            Dg = np.insert(
                atm.Dzz[:, species.index('H2SO4')], 0, atm.Dzz[0, species.index('H2SO4')]
            )
            rate = (
                Dg
                * m
                / rho_p
                / r_p**2
                * (var.y[:, species.index('H2SO4')] - atm.sat_p['H2SO4'] / kb / atm.Tco)
            )

            # rate_c = m/(4*rho_p)*(8*kb*atm.Tco/np.pi/m)**0.5 *(var.y[:,species.index('H2SO4')]-atm.sat_p['H2SO4']/kb/atm.Tco)/r_p

            # how many gas molecules are contained in one particle with the assumed size r_p
            n_mol = 4.0 / 3 * np.pi * r_p**3 * rho_p / m

            var.k[re] = rate
            var.k[re + 1] = rate  # /n_mol

            # positive: condensation
            var.k[re] = np.maximum(var.k[re], 0)
            # negative: evaporation
            var.k[re + 1] = np.minimum(var.k[re + 1], 0)
            var.k[re + 1] = np.abs(var.k[re + 1])

        elif var.Rf[re] == 'S2 -> S2_l_s' and 'S2' in self.cfg.condense_sp:
            m = 45.019 / Navo
            rho_p = atm.rho_p['S2_l_s']
            r_p = atm.r_p['S2_l_s']

            # new approach: contiuum regime DM/rho c
            Dg = np.insert(
                atm.Dzz[:, species.index('S2')], 0, atm.Dzz[0, species.index('S2')]
            )
            rate = (
                Dg
                * m
                / rho_p
                / r_p**2
                * (var.y[:, species.index('S2')] - atm.sat_p['S2'] / kb / atm.Tco)
            )

            # how many gas molecules are contained in one particle with the assumed size r_p
            n_mol = 4.0 / 3 * np.pi * r_p**3 * rho_p / m

            var.k[re] = rate
            var.k[re + 1] = rate  # /n_mol

            # positive: condensation
            var.k[re] = np.maximum(var.k[re], 0)
            # negative: evaporation
            var.k[re + 1] = np.minimum(var.k[re + 1], 0)
            var.k[re + 1] = np.abs(var.k[re + 1])

        elif var.Rf[re] == 'S4 -> S4_l_s' and 'S4' in self.cfg.condense_sp:
            m = 32.06 * 4 / Navo
            rho_p = atm.rho_p['S4_l_s']
            r_p = atm.r_p['S4_l_s']

            # new approach: contiuum regime DM/rho c
            Dg = np.insert(
                atm.Dzz[:, species.index('S4')], 0, atm.Dzz[0, species.index('S4')]
            )
            rate = (
                Dg
                * m
                / rho_p
                / r_p**2
                * (var.y[:, species.index('S4')] - atm.sat_p['S4'] / kb / atm.Tco)
            )

            # how many gas molecules are contained in one particle with the assumed size r_p
            n_mol = 4.0 / 3 * np.pi * r_p**3 * rho_p / m

            # acc_ratio = var.y[:,species.index('S4_l_s')]/n_mol /self.cfg.n_ccn # accomdation ratio: 0 all ccn available 1=no more free ccn
            # lim_factor = 1-acc_ratio
            # lim_factor[lim_factor<0] = 0

            var.k[re] = rate  # *lim_factor
            var.k[re + 1] = rate  # /n_mol

            # positive: condensation
            var.k[re] = np.maximum(var.k[re], 0)
            # negative: evaporation
            var.k[re + 1] = np.minimum(var.k[re + 1], 0)
            var.k[re + 1] = np.abs(var.k[re + 1])

        elif var.Rf[re] == 'S8 -> S8_l_s' and 'S8' in self.cfg.condense_sp:
            m = 360.152 / Navo
            rho_p = atm.rho_p['S8_l_s']
            r_p = atm.r_p['S8_l_s']

            # new approach: contiuum regime DM/rho c
            Dg = np.insert(
                atm.Dzz[:, species.index('S8')], 0, atm.Dzz[0, species.index('S8')]
            )
            rate = (
                Dg
                * m
                / rho_p
                / r_p**2
                * (var.y[:, species.index('S8')] - atm.sat_p['S8'] / kb / atm.Tco)
            )

            # how many gas molecules are contained in one particle with the assumed size r_p
            n_mol = 4.0 / 3 * np.pi * r_p**3 * rho_p / m

            var.k[re] = rate
            var.k[re + 1] = rate  # /n_mol

            # positive: condensation
            var.k[re] = np.maximum(var.k[re], 0)
            # negative: evaporation
            var.k[re + 1] = np.minimum(var.k[re + 1], 0)
            var.k[re + 1] = np.abs(var.k[re + 1])

        elif var.Rf[re] == 'C -> C_s' and 'C' in self.cfg.condense_sp:
            m = 12.011 / Navo
            rho_p = atm.rho_p['C_s']
            r_p = atm.r_p['C_s']

            # new approach: contiuum regime DM/rho c
            Dg = np.insert(
                atm.Dzz[:, species.index('C')], 0, atm.Dzz[0, species.index('C')]
            )
            rate = (
                Dg
                * m
                / rho_p
                / r_p**2
                * (var.y[:, species.index('C')] - atm.sat_p['C'] / kb / atm.Tco)
            )

            # how many gas molecules are contained in one particle with the assumed size r_p
            n_mol = 4.0 / 3 * np.pi * r_p**3 * rho_p / m

            var.k[re] = rate
            var.k[re + 1] = rate  # /n_mol

            # positive: condensation
            var.k[re] = np.maximum(var.k[re], 0)
            # negative: evaporation
            var.k[re + 1] = np.minimum(var.k[re + 1], 0)
            var.k[re + 1] = np.abs(var.k[re + 1])

    return var

`conv(var, para, atm, out=False, print_freq=100)`

funtion returns TRUE if the convergence condition is satisfied

Source code in src/vulcan/op.py

def conv(self, var, para, atm, out=False, print_freq=100):
    """
    funtion returns TRUE if the convergence condition is satisfied
    """
    st_factor, mtol_conv, atol, yconv_cri, slope_cri, yconv_min = (
        self.cfg.st_factor,
        self.cfg.mtol_conv,
        self.cfg.atol,
        self.cfg.yconv_cri,
        self.cfg.slope_cri,
        self.cfg.yconv_min,
    )
    y, ymix, y_time, t_time = var.y.copy(), var.ymix.copy(), var.y_time, var.t_time
    count = para.count

    # slope_min = min( np.amin(atm.Kzz)/np.amax(0.1*atm.Hp)**2 , 1.e-8)
    slope_min = min(np.amin(atm.Kzz / (0.1 * atm.Hp[:-1]) ** 2), 1.0e-8)
    slope_min = max(slope_min, 1.0e-10)

    indx = np.abs(t_time - var.t * st_factor).argmin()
    if indx == para.count - 1:
        indx -= 1  # Important!! For dt larger than half of the runtime (count-1 is the last one)

    # Don't check more than self.cfg.conv_step (1000) steps back
    indx = max(para.count - self.cfg.conv_step, indx)

    # TEST
    longdy = np.abs((y_time[count - 1] - y_time[indx]) / np.vstack(atm.n_0))
    longdy[ymix < mtol_conv] = 0
    longdy[y < atol] = 0

    # to get rid off non-convergent species, e.g. HC3N without sinks
    if 'conver_ignore' in dir(self.cfg):
        for sp in self.cfg.conver_ignore:
            longdy[:, species.index(sp)] = 0  # added 2023

    if self.cfg.use_condense:
        longdy[:, self.non_gas_sp_index] = 0

    with np.errstate(
        divide='ignore', invalid='ignore'
    ):  # ignoring nan when devided by zero
        where_varies_most = longdy / ymix
    para.where_varies_most = where_varies_most

    longdy = np.amax(longdy[ymix > 0] / ymix[ymix > 0])
    longdydt = longdy / (t_time[-1] - t_time[indx])
    # store longdy and longdydt
    var.longdy, var.longdydt = longdy, longdydt

    if (
        longdy < yconv_cri
        and longdydt < slope_cri
        or longdy < yconv_min
        and longdydt < slope_min
    ) and var.aflux_change < self.cfg.flux_cri:
        return True

    return False

`f_dy(var, para)`

Source code in src/vulcan/op.py

def f_dy(self, var, para):  # y, y_prev, ymix, yconv, count, dt
    if para.count == 0:
        var.dy, var.dydt = 1.0, 1.0
        return var
    y, ymix, y_prev = var.y, var.ymix, var.y_prev
    dy = np.abs(y - y_prev)
    dy[ymix < self.cfg.mtol] = 0
    dy[y < self.cfg.atol] = 0
    dy = np.amax(dy[y > 0] / y[y > 0])

    var.dy, var.dydt = dy, dy / var.dt

    return var

`h2o_conden_evap_relax(var, atm)`

Source code in src/vulcan/op.py

def h2o_conden_evap_relax(self, var, atm):
    m = 18.0 / Navo
    rho_p = atm.rho_p['H2O_l_s']  # mix of water and ice
    r_p = atm.r_p['H2O_l_s']
    # relative humidity
    sat_humidity = atm.sat_p['H2O'] / kb / atm.Tco * self.cfg.humidity

    # new approach: contiuum regime DM/rho c
    Dg = np.insert(atm.Dzz[:, species.index('H2O')], 0, atm.Dzz[0, species.index('H2O')])
    tau = 1.0 / (Dg * m / rho_p / r_p**2 * (var.y[:, species.index('H2O')] - sat_humidity))
    conden_indx = np.where(tau > 0)
    evap_indx = np.where(tau < 0)
    sat_mix = sat_humidity / atm.n_0
    # tau = np.abs(tau)

    # implicit-Euler to remove water
    y_conden = (var.ymix[:, species.index('H2O')] + var.dt / tau * sat_mix) / (
        1.0 + var.dt / tau
    )

    # evaporation to remove ice/water(liquid)
    ice_loss = (
        (var.y[:, species.index('H2O')] - sat_humidity) * var.dt / tau
    )  # both tau < 0 and y_H2O - sat < 0
    # cannot lose more than it has
    ice_loss = np.minimum(var.y[:, species.index('H2O_l_s')], ice_loss)

    # how many gas molecules are contained in one particle with the assumed size r_p
    # n_mol = 4./3*np.pi*r_p**3 *rho_p /m
    # and converting the mixing ratio of molecules /cm3 to droplets/cm3
    # "move" the condensed water to H2O_l_s
    var.ymix[conden_indx, species.index('H2O_l_s')] += (
        var.ymix[conden_indx, species.index('H2O')] - y_conden[conden_indx]
    )
    var.ymix[conden_indx, species.index('H2O')] = y_conden[conden_indx]
    # store the saturated parts (only relax where ymix > ysat)

    var.ymix[evap_indx, species.index('H2O')] += ice_loss[evap_indx] / atm.n_0[evap_indx]
    var.ymix[evap_indx, species.index('H2O_l_s')] -= (
        ice_loss[evap_indx] / atm.n_0[evap_indx]
    )

    var.y = var.ymix * np.vstack(np.sum(var.y[:, atm.gas_indx], axis=1))

    return var

`h2o_conden_relax(var, atm)`

Source code in src/vulcan/op.py

def h2o_conden_relax(self, var, atm):
    m = 18.0 / Navo
    rho_p = 0.95  # mix of water and ice
    r_p = atm.r_p['H2O_l_s']
    # relative humidity
    sat_humidity = atm.sat_p['H2O'] / kb / atm.Tco * self.cfg.humidity

    # new approach: contiuum regime DM/rho c
    Dg = np.insert(atm.Dzz[:, species.index('H2O')], 0, atm.Dzz[0, species.index('H2O')])
    tau = np.abs(
        1.0 / (Dg * m / rho_p / r_p**2 * (var.y[:, species.index('H2O')] - sat_humidity))
    )
    sat_mix = sat_humidity / atm.n_0

    # implicit-Euler to remove water
    y_conden = (var.ymix[:, species.index('H2O')] + var.dt / tau * sat_mix) / (
        1.0 + var.dt / tau
    )
    conden_indx = np.where(var.ymix[:, species.index('H2O')] > sat_mix)

    # how many gas molecules are contained in one particle with the assumed size r_p
    n_mol = 4.0 / 3 * np.pi * r_p**3 * rho_p / m
    # and converting the mixing ratio of molecules /cm3 to droplets/cm3
    # "move" the condensed water to H2O_l_s
    var.ymix[conden_indx, species.index('H2O_l_s')] += (
        var.ymix[conden_indx, species.index('H2O')] - y_conden[conden_indx]
    )  # /n_mol

    var.ymix[conden_indx, species.index('H2O')] = y_conden[conden_indx]
    # restore the unsaturated parts (only relax where ymix > ysat)

    var.y = var.ymix * np.vstack(np.sum(var.y[:, atm.gas_indx], axis=1))

    return var

`nh3_conden_evap_relax(var, atm)`

Source code in src/vulcan/op.py

def nh3_conden_evap_relax(self, var, atm):
    m = 17.0 / Navo
    rho_p = atm.rho_p['NH3_l_s']  # mix of water and ice
    r_p = atm.r_p['NH3_l_s']
    # relative humidity
    sat_p = atm.sat_p['NH3'] / kb / atm.Tco
    sat_mix = sat_p / atm.n_0

    conden_top = np.argmin(sat_mix)

    # new approach: contiuum regime DM/rho c
    Dg = np.insert(atm.Dzz[:, species.index('NH3')], 0, atm.Dzz[0, species.index('NH3')])
    tau = 1.0 / (Dg * m / rho_p / r_p**2 * (var.y[:, species.index('NH3')] - sat_p))
    conden_indx = np.where(tau > 0)[0]
    evap_indx = np.where(tau < 0)[0]

    # above the top of condensation zone, there should NOT be any condensation when using the relaxiation method
    conden_indx = [i for i in conden_indx if i <= conden_top]
    # evap_indx = [i for i in evap_indx if i <= conden_top]

    # implicit-Euler to remove water
    y_conden = (var.ymix[:, species.index('NH3')] + var.dt / tau * sat_mix) / (
        1.0 + var.dt / tau
    )

    # evaporation to remove ice/water(liquid)
    ice_loss = (
        (var.y[:, species.index('NH3')] - sat_p) * var.dt / tau
    )  # both tau < 0 and y_H2O - sat < 0
    # cannot lose more than it has
    ice_loss = np.minimum(var.y[:, species.index('NH3_l_s')], ice_loss)

    # how many gas molecules are contained in one particle with the assumed size r_p
    # n_mol = 4./3*np.pi*r_p**3 *rho_p /m
    # and converting the mixing ratio of molecules /cm3 to droplets/cm3
    # "move" the condensed water to H2O_l_s
    var.ymix[conden_indx, species.index('NH3_l_s')] += (
        var.ymix[conden_indx, species.index('NH3')] - y_conden[conden_indx]
    )
    var.ymix[conden_indx, species.index('NH3')] = y_conden[conden_indx]
    # store the saturated parts (only relax where ymix > ysat)

    var.ymix[evap_indx, species.index('NH3')] += ice_loss[evap_indx] / atm.n_0[evap_indx]
    # instaneous evaporation
    var.ymix[evap_indx, species.index('NH3_l_s')] -= (
        ice_loss[evap_indx] / atm.n_0[evap_indx]
    )
    # var.ymix[evap_indx,species.index('NH3_l_s')] = 0

    var.ymix[:, species.index('NH3_l_s')] = np.maximum(
        var.ymix[:, species.index('NH3_l_s')], 0
    )  # cannot lose more than it has

    var.y = var.ymix * np.vstack(np.sum(var.y[:, atm.gas_indx], axis=1))

    return var

`save_step(var, para)`

save current values of y and add 1 to the counter

Source code in src/vulcan/op.py

def save_step(self, var, para):
    """
    save current values of y and add 1 to the counter
    """
    var.t += var.dt
    para.count += 1

    # tmp = list(var.y)
    # if para.count % self.y_time_freq ==0:
    var.y_time.append(var.y)
    # var.ymix_time.append(var.ymix.copy())
    var.t_time.append(var.t)

    # only used in PI_control
    # var.dy_time.append(var.y)
    # var.dydt_time.append(var.dydt)
    var.atom_loss_time.append(list(var.atom_loss.values()))

    return var, para

`stop(var, para, atm)`

To check the convergence criteria and stop the integration

Source code in src/vulcan/op.py

def stop(self, var, para, atm):
    """
    To check the convergence criteria and stop the integration
    """
    if (
        var.t > self.cfg.trun_min
        and para.count > self.cfg.count_min
        and self.conv(var, para, atm)
    ):
        log.info('Integration successful with ' + str(para.count) + ' steps')
        log.info('long dy, long dydt = ' + str(var.longdy) + ', ' + str(var.longdydt))
        log.debug('Actinic flux change: ' + '{:.2E}'.format(var.aflux_change))
        self.output.print_end_msg(var, para)
        para.end_case = 1
        return True
    elif var.t > self.cfg.runtime:
        log.warning(
            'After ------- %s seconds -------' % (time.time() - para.start_time)
            + ' s CPU time'
        )
        log.warning('Integration not completed...')
        log.warning('Maximal allowed runtime exceeded (' + str(self.cfg.runtime) + ' sec)!')
        para.end_case = 2
        return True
    elif para.count > self.cfg.count_max:
        log.warning(
            'After ------- %s seconds -------' % (time.time() - para.start_time)
            + ' s CPU time'
        )
        log.warning('Integration not completed...')
        log.warning('Maximal allowed steps exceeded (' + str(self.cfg.count_max) + ')!')
        para.end_case = 3
        return True

`update_mu_dz(var, atm, make_atm)`

Source code in src/vulcan/op.py

def update_mu_dz(self, var, atm, make_atm):  # y, ni, spec, Tco, pco

    # gravity
    gz = atm.g
    pref_indx = atm.pref_indx
    Tco, pico = atm.Tco, atm.pico.copy()
    # calculating mu (mean molecular weight)
    atm = make_atm.mean_mass(var, atm, ni)
    Hp = atm.Hp
    nz = self.cfg.nz

    for i in range(pref_indx, nz):
        if i == pref_indx:
            atm.g[i] = atm.gs
            Hp[i] = kb * Tco[i] / (atm.mu[i] / Navo * atm.gs)
        else:
            atm.g[i] = atm.gs * (self.cfg.Rp / (self.cfg.Rp + atm.zco[i])) ** 2
            Hp[i] = kb * Tco[i] / (atm.mu[i] / Navo * atm.g[i])
        atm.dz[i] = Hp[i] * np.log(
            pico[i] / pico[i + 1]
        )  # pico[i+1] has a lower P than pico[i] (higer height)
        atm.zco[i + 1] = atm.zco[i] + atm.dz[i]  # zco is set zero at 1bar for gas giants

    # for pref_indx != zero
    if not pref_indx == 0:
        for i in range(pref_indx - 1, -1, -1):
            atm.g[i] = atm.gs * (self.cfg.Rp / (self.cfg.Rp + atm.zco[i + 1])) ** 2
            Hp[i] = kb * Tco[i] / (atm.mu[i] / Navo * atm.g[i])
            atm.dz[i] = Hp[i] * np.log(pico[i] / pico[i + 1])
            atm.zco[i] = atm.zco[i + 1] - atm.dz[i]  # from i+1 propogating down to i

    zmco = 0.5 * (atm.zco + np.roll(atm.zco, -1))
    atm.zmco = zmco[:-1]
    dzi = 0.5 * (atm.dz + np.roll(atm.dz, 1))
    atm.dzi = dzi[1:]

    # for the molecular diffsuion
    if self.cfg.use_moldiff:
        Ti = 0.5 * (Tco + np.roll(Tco, -1))
        atm.Ti = Ti[:-1]
        Hpi = 0.5 * (Hp + np.roll(Hp, -1))
        atm.Hpi = Hpi[:-1]

    return atm

`update_phi_esc(var, atm)`

Source code in src/vulcan/op.py

def update_phi_esc(self, var, atm):  # updating diffusion-mimited escape

    # Diffusion limited escape
    for sp in self.cfg.diff_esc:
        # atm.top_flux[species.index(sp)] = - atm.Dzz[-1,species.index(sp)] *var.y[-1,species.index(sp)] /atm.Hp[-1]
        atm.top_flux[species.index(sp)] = (
            -atm.Dzz[-1, species.index(sp)]
            * var.y[-1, species.index(sp)]
            * (
                1.0 / atm.Hp[-1]
                - atm.ms[species.index(sp)] * atm.g[-1] / (Navo * kb * atm.Tco[-1])
            )
        )
        atm.top_flux[species.index(sp)] = max(
            atm.top_flux[species.index(sp)], self.cfg.max_flux * (-1)
        )

    return atm

`ODESolver(vulcan_cfg)`

Bases: object

Source code in src/vulcan/op.py

def __init__(self, vulcan_cfg: Config):  # do I always need to update var, atm, para ?

    self.cfg = vulcan_cfg
    self.mtol = self.cfg.mtol
    self.atol = self.cfg.atol
    self.non_gas_sp = self.cfg.non_gas_sp

    if self.cfg.use_condense:
        self.non_gas_sp_index = [species.index(sp) for sp in self.non_gas_sp]
        self.condense_sp_index = [species.index(sp) for sp in self.cfg.condense_sp]

    self.fix_sp_bot_index = [species.index(sp) for sp in self.cfg.use_fix_sp_bot.keys()]
    self.fix_sp_bot_mix = np.array(
        [self.cfg.use_fix_sp_bot[sp] for sp in self.cfg.use_fix_sp_bot.keys()]
    )

`atol = self.cfg.atol` `instance-attribute`

`cfg = vulcan_cfg` `instance-attribute`

`condense_sp_index = [(species.index(sp)) for sp in (self.cfg.condense_sp)]` `instance-attribute`

`fix_sp_bot_index = [(species.index(sp)) for sp in (self.cfg.use_fix_sp_bot.keys())]` `instance-attribute`

`fix_sp_bot_mix = np.array([(self.cfg.use_fix_sp_bot[sp]) for sp in (self.cfg.use_fix_sp_bot.keys())])` `instance-attribute`

`mtol = self.cfg.mtol` `instance-attribute`

`non_gas_sp = self.cfg.non_gas_sp` `instance-attribute`

`non_gas_sp_index = [(species.index(sp)) for sp in (self.non_gas_sp)]` `instance-attribute`

`clip(var, para, atm, pos_cut=0, nega_cut=-1)`

function to clip samll and negative values and to calculate the particle loss

Source code in src/vulcan/op.py

def clip(self, var, para, atm, pos_cut=0, nega_cut=-1):
    """
    function to clip samll and negative values
    and to calculate the particle loss
    """
    y, ymix = var.y, var.ymix.copy()

    para.small_y += np.abs(np.sum(y[np.logical_and(y < pos_cut, y >= 0)]))
    para.nega_y += np.abs(np.sum(y[np.logical_and(y > nega_cut, y <= 0)]))
    y[np.logical_and(y < pos_cut, y >= nega_cut)] = 0.0

    # Also setting y=0 when ymix<mtol
    y[np.logical_and(ymix < self.mtol, y < 0)] = 0.0

    var = self.loss(var)

    # store y and ymix
    # TEST condensation excluding non-gaseous species
    if self.cfg.non_gas_sp:
        var.y, var.ymix = y, var.y / np.vstack(np.sum(var.y[:, atm.gas_indx], axis=1))
    else:
        var.y, var.ymix = y, y / np.vstack(np.sum(y, axis=1))
    # TEST condensation excluding non-gaseous species

    return var, para

`compute_J(var, atm)`

computes photodissociation/photoionization rates; including T-dependent cross sections

Source code in src/vulcan/op.py

def compute_J(self, var, atm):  # the vectorized version
    """
    computes photodissociation/photoionization rates; including T-dependent cross sections
    """
    flux = var.aflux
    nz = self.cfg.nz

    diss_cross = var.cross_J  # use the key (sp, branch index) e.g. ("H2O", 1); 1D array
    diss_cross_T = var.cross_J_T  # 2D array with the shape of nz * bins

    bins = var.bins
    n_branch = var.n_branch

    # reset to zeros every time
    var.J_sp = dict(
        [((sp, bn), np.zeros(nz)) for sp in var.photo_sp for bn in range(n_branch[sp] + 1)]
    )

    for sp in var.photo_sp:
        # shape: flux (nz,nbin) cross (nbin)

        for nbr in range(1, n_branch[sp] + 1):  # axis=1 is to sum over all wavelength
            if sp in self.cfg.T_cross_sp:
                var.J_sp[(sp, nbr)] = np.sum(
                    flux[:, : var.sflux_din12_indx]
                    * diss_cross_T[(sp, nbr)][:, : var.sflux_din12_indx]
                    * var.dbin1,
                    axis=1,
                )
                var.J_sp[(sp, nbr)] -= (
                    0.5
                    * (
                        flux[:, 0] * diss_cross_T[(sp, nbr)][:, 0]
                        + flux[:, var.sflux_din12_indx - 1]
                        * diss_cross_T[(sp, nbr)][:, var.sflux_din12_indx - 1]
                    )
                    * var.dbin1
                )
                var.J_sp[(sp, nbr)] += np.sum(
                    flux[:, var.sflux_din12_indx :]
                    * diss_cross_T[(sp, nbr)][:, var.sflux_din12_indx :]
                    * var.dbin2,
                    axis=1,
                )
                var.J_sp[(sp, nbr)] -= (
                    0.5
                    * (
                        flux[:, var.sflux_din12_indx]
                        * diss_cross_T[(sp, nbr)][:, var.sflux_din12_indx]
                        + flux[:, -1] * diss_cross_T[(sp, nbr)][:, -1]
                    )
                    * var.dbin2
                )

            else:
                var.J_sp[(sp, nbr)] = np.sum(
                    flux[:, : var.sflux_din12_indx]
                    * diss_cross[(sp, nbr)][: var.sflux_din12_indx]
                    * var.dbin1,
                    axis=1,
                )
                var.J_sp[(sp, nbr)] -= (
                    0.5
                    * (
                        flux[:, 0] * diss_cross[(sp, nbr)][0]
                        + flux[:, var.sflux_din12_indx - 1]
                        * diss_cross[(sp, nbr)][var.sflux_din12_indx - 1]
                    )
                    * var.dbin1
                )
                var.J_sp[(sp, nbr)] += np.sum(
                    flux[:, var.sflux_din12_indx :]
                    * diss_cross[(sp, nbr)][var.sflux_din12_indx :]
                    * var.dbin2,
                    axis=1,
                )
                var.J_sp[(sp, nbr)] -= (
                    0.5
                    * (
                        flux[:, var.sflux_din12_indx]
                        * diss_cross[(sp, nbr)][var.sflux_din12_indx]
                        + flux[:, -1] * diss_cross[(sp, nbr)][-1]
                    )
                    * var.dbin2
                )

            # summing over all branches
            var.J_sp[(sp, 0)] += var.J_sp[(sp, nbr)]
            # incoperating J into rate coefficients
            if var.pho_rate_index[(sp, nbr)] not in self.cfg.remove_list:
                var.k[var.pho_rate_index[(sp, nbr)]] = (
                    var.J_sp[(sp, nbr)] * self.cfg.f_diurnal
                )  # f_diurnal = 0.5 for Earth; = 1 for tidally-loced planets

`compute_Jion(var, atm)`

compute the photoionization rate haven't considered any temperature dependence yet

Source code in src/vulcan/op.py

def compute_Jion(self, var, atm):
    """
    compute the photoionization rate
    haven't considered any temperature dependence yet
    """
    nz = self.cfg.nz
    flux = var.aflux
    ion_cross = var.cross_Jion  # use the key (sp, br) e.g. ("H2O", 1)

    bins = var.bins
    n_branch = var.ion_branch

    # reset to zeros every time
    var.Jion_sp = dict(
        [((sp, bn), np.zeros(nz)) for sp in var.ion_sp for bn in range(n_branch[sp] + 1)]
    )

    for sp in var.ion_sp:
        # shape: flux (nz,nbin) cross (nbin)

        # convert to actinic flux *1/(hc/ld)
        for nbr in range(1, n_branch[sp] + 1):
            # axis=1 is to sum over all wavelength
            var.Jion_sp[(sp, nbr)] = np.sum(
                flux[:, : var.sflux_din12_indx]
                * ion_cross[(sp, nbr)][: var.sflux_din12_indx]
                * var.dbin1,
                axis=1,
            )
            var.Jion_sp[(sp, nbr)] -= (
                0.5
                * (
                    flux[:, 0] * ion_cross[(sp, nbr)][0]
                    + flux[:, var.sflux_din12_indx - 1]
                    * ion_cross[(sp, nbr)][var.sflux_din12_indx - 1]
                )
                * var.dbin1
            )
            var.Jion_sp[(sp, nbr)] += np.sum(
                flux[:, var.sflux_din12_indx :]
                * ion_cross[(sp, nbr)][var.sflux_din12_indx :]
                * var.dbin2,
                axis=1,
            )
            var.Jion_sp[(sp, nbr)] -= (
                0.5
                * (
                    flux[:, var.sflux_din12_indx]
                    * ion_cross[(sp, nbr)][var.sflux_din12_indx]
                    + flux[:, -1] * ion_cross[(sp, nbr)][-1]
                )
                * var.dbin2
            )

            # 0 is the total dissociation rate
            # summing all branches

            var.Jion_sp[(sp, 0)] += var.Jion_sp[(sp, nbr)]
            # incoperating J into rate coefficients
            if var.ion_rate_index[(sp, nbr)] not in self.cfg.remove_list:
                var.k[var.ion_rate_index[(sp, nbr)]] = (
                    var.Jion_sp[(sp, nbr)] * self.cfg.f_diurnal
                )  # f_diurnal = 0.5 for Earth; = 1 for tidally-loced planets

`compute_flux(var, atm)`

Source code in src/vulcan/op.py

def compute_flux(self, var, atm):  # Vectorise this loop!
    # change it to stagerred grids
    # top: stellar flux
    # bottom BC: zero upcoming flux

    # Note!!! Matej's mu is defined in the outgoing hemisphere so his mu<0
    # My cos[sl_angle] is always 0<=mu<=1
    # Converting my mu to Matej's mu (e.g. 45 deg -> 135 deg)

    nz = self.cfg.nz

    mu_ang = -1.0 * np.cos(self.cfg.sl_angle)
    edd = self.cfg.edd
    tau = var.tau

    # delta_tau (length nz) is used in the transmission function
    delta_tau = tau - np.roll(
        tau, -1, axis=0
    )  # np.roll(tau,-1,axis=0) are the upper layers
    delta_tau = delta_tau[:-1]

    # single-scattering albedo
    nbins = len(var.bins)
    tot_abs, tot_scat = np.zeros((nz, nbins)), np.zeros((nz, nbins))
    for sp in var.photo_sp:
        tot_abs += np.vstack(var.ymix[:, species.index(sp)]) * var.cross[sp]  # nz * nbins
    for sp in self.cfg.scat_sp:
        tot_scat += np.vstack(var.ymix[:, species.index(sp)]) * var.cross_scat[sp]

    total = tot_abs + tot_scat

    w0 = tot_scat / (tot_abs + tot_scat)  # 2D: nz * nbins
    # tot_abs + tot_scat can be zero when certain gas (e.g. H2) does not exist

    # Replace nan with zero and inf with very large numbers
    w0 = np.nan_to_num(w0)

    # to avoit w0=1
    w0 = np.minimum(w0, 1.0 - 1.0e-8)

    # sflux: the direct beam; dflux: diffusive flux
    """ Beer's law for the intensity"""
    var.sflux = var.sflux_top * np.exp(-1.0 * tau / np.cos(self.cfg.sl_angle))
    # converting the intensity to flux for the raditive transfer calculation
    dir_flux = (
        var.sflux * np.cos(self.cfg.sl_angle)
    )  # need to convert to diffuse flux in the RT definition so it can covert back to total intensity with eps

    # scattering
    # the transmission function (length nz)
    if ag0 == 0:  # to save memory
        tran = np.exp(-1.0 / edd * (1.0 - w0) ** 0.5 * delta_tau)  # 2D: nz * nbins
        zeta_p = 0.5 * (1.0 + (1.0 - w0) ** 0.5)
        zeta_m = 0.5 * (1.0 - (1.0 - w0) ** 0.5)
        ll = -1.0 * w0 / (1.0 / mu_ang**2 - 1.0 / edd**2 * (1.0 - w0))
        g_p = 0.5 * (ll * (1.0 / edd + 1.0 / mu_ang))
        g_m = 0.5 * (ll * (1.0 / edd - 1.0 / mu_ang))

    else:
        tran = np.exp(-1.0 / edd * ((1.0 - w0 * ag0) * (1.0 - w0)) ** 0.5 * delta_tau)
        zeta_p = 0.5 * (1.0 + ((1.0 - w0) / (1 - w0 * ag0)) ** 0.5)
        zeta_m = 0.5 * (1.0 - ((1.0 - w0) / (1 - w0 * ag0)) ** 0.5)
        ll = ((1.0 - w0) * (1 - w0 * ag0) - 1.0) / (
            1.0 / mu_ang**2 - 1.0 / edd**2 * (1.0 - w0) * (1 - w0 * ag0)
        )
        g_p = 0.5 * (
            ll * (1.0 / edd + 1 / (mu_ang * (1.0 - w0 * ag0)))
            + w0 * ag0 * mu_ang / (1.0 - w0 * ag0)
        )
        g_m = 0.5 * (
            ll * (1.0 / edd - 1 / (mu_ang * (1.0 - w0 * ag0)))
            - w0 * ag0 * mu_ang / (1.0 - w0 * ag0)
        )

    # to avoit zero denominator
    ll = np.minimum(ll, 1.0e10)
    ll = np.maximum(ll, -1.0e10)

    # 2D: nz * nbins
    chi = zeta_m**2 * tran**2 - zeta_p**2
    xi = zeta_p * zeta_m * (1.0 - tran**2)
    phi = (zeta_m**2 - zeta_p**2) * tran

    # 2D: nz * nbins
    i_u = phi * g_p * dir_flux[:-1] - (xi * g_m + chi * g_p) * dir_flux[1:]
    i_d = phi * g_m * dir_flux[1:] - (chi * g_m + xi * g_p) * dir_flux[:-1]
    # sflux[1:] are all the layers above and sflux[:-1] are all the layers abelow

    var.zeta_m = zeta_m
    var.zeta_p = zeta_p
    var.tran = tran

    # For testing computating speed
    # starting recording time
    # start_time = timeit.default_timer()

    # propagating downward layer by layer and then upward
    # var.dflux_d and var.dflux_p are defined at the interfaces (staggerred)
    # the rest is defined in the center of the layer
    for j in range(
        nz - 1, -1, -1
    ):  # dflux_d goes from the second top interface (nz+1 interfaces)
        var.dflux_d[j] = (
            1.0
            / chi[j]
            * (phi[j] * var.dflux_d[j + 1] - xi[j] * var.dflux_u[j] + i_d[j] / mu_ang)
        )
    for j in range(1, nz + 1):
        var.dflux_u[j] = (
            1.0
            / chi[j - 1]
            * (
                phi[j - 1] * var.dflux_u[j - 1]
                - xi[j - 1] * var.dflux_d[j]
                + i_u[j - 1] / mu_ang
            )
        )

    # the average flux from the direct beam
    # !!! WITHOUT multiplied by the cos zenith angle (flux per unit area perpendicular to the direction of propagationat) !!!
    ave_dir_flux = 0.5 * (var.sflux[:-1] + var.sflux[1:])
    # devided by the Eddington coefficient to recover the total intensity (integrated over all directions)
    tot_flux = (
        ave_dir_flux
        + 0.5
        * (var.dflux_u[:-1] + var.dflux_u[1:] + var.dflux_d[1:] + var.dflux_d[:-1])
        / edd
    )

    # store the previous actinic flux into prev_aflux
    var.prev_aflux = np.copy(var.aflux)
    # converting to the actinic flux and storing the current flux
    var.aflux = tot_flux / (hc / var.bins)
    # the change of the actinic flux
    var.aflux_change = np.nanmax(
        np.abs(var.aflux - var.prev_aflux)[var.aflux > self.cfg.flux_atol]
        / var.aflux[var.aflux > self.cfg.flux_atol]
    )

`compute_tau(var, atm)`

compute the optical depth

Source code in src/vulcan/op.py

def compute_tau(self, var, atm):
    """compute the optical depth"""

    nz = self.cfg.nz

    # reset to zero
    var.tau.fill(0)
    # absorption species
    absp_sp = set.union(var.photo_sp, var.ion_sp)

    for j in range(nz - 1, -1, -1):
        for sp in absp_sp:
            # summing over all T-dependentphoto species
            if sp in self.cfg.T_cross_sp:
                var.tau[j] += (
                    var.y[j, species.index(sp)] * atm.dz[j] * var.cross_T[sp][j]
                )  # 1-D shape of nbins from the j level
            else:  # summing over all T-independent photo species
                var.tau[j] += (
                    var.y[j, species.index(sp)] * atm.dz[j] * var.cross[sp]
                )  # only the j-th laye

        for sp in self.cfg.scat_sp:  # scat_sp are not necessary photo_sp, e.g. He
            var.tau[j] += var.y[j, species.index(sp)] * atm.dz[j] * var.cross_scat[sp]
        # adding the layer above at the end of species loop
        var.tau[j] += var.tau[j + 1]

`diffdf(y, atm)`

function of eddy diffusion including molecular diffusion, with zero-flux boundary conditions and non-uniform grids (dzi) in the form of Ajy_j + Bj+1y_j+1 + Cj-1*y_j-1

Source code in src/vulcan/op.py

def diffdf(self, y, atm):
    """
    function of eddy diffusion including molecular diffusion, with zero-flux boundary conditions and non-uniform grids (dzi)
    in the form of Aj*y_j + Bj+1*y_j+1 + Cj-1*y_j-1
    """

    nz = self.cfg.nz

    y = y.copy()

    # TEST condensation excluding non-gaseous species
    if self.cfg.non_gas_sp:
        ysum = np.sum(y[:, atm.gas_indx], axis=1)
    else:
        ysum = np.sum(y, axis=1)
    # TEST condensation excluding non-gaseous species

    dzi = atm.dzi.copy()
    Kzz = atm.Kzz.copy()
    vz = atm.vz.copy()
    Dzz = atm.Dzz.copy()
    alpha = atm.alpha.copy()
    Tco = atm.Tco
    ms = atm.ms.copy()
    Hp = atm.Hp.copy()
    g = atm.g
    Ti = atm.Ti
    Hpi = atm.Hpi

    # # define T_1/2 for the molecular diffusion
    # Ti = 0.5*(Tco + np.roll(Tco,-1))
    # Ti = Ti[:-1]
    # Hpi = 0.5*(Hp + np.roll(Hp,-1))
    # Hpi = Hpi[:-1]
    # # store Ti and Hpi
    # atm.Ti = Ti
    # atm.Hpi = Hpi

    A, B, C = np.zeros(nz), np.zeros(nz), np.zeros(nz)
    Ai, Bi, Ci = [np.zeros((nz, ni)) for i in range(3)]

    A[0] = -1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / 2.0 / ysum[0]
    B[0] = 1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / 2.0 / ysum[1]
    C[0] = 0
    A[nz - 1] = (
        -1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[nz - 1] + ysum[nz - 2])
        / 2.0
        / ysum[nz - 1]
    )
    B[nz - 1] = 0
    C[nz - 1] = (
        1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[nz - 1] + ysum[nz - 2])
        / 2.0
        / ysum[nz - 2]
    )

    # vertical adection (with closed B.C.)
    A[0] += -((vz[0] > 0) * vz[0]) / dzi[0]
    B[0] += -((vz[0] < 0) * vz[0]) / dzi[0]
    A[-1] += ((vz[-1] < 0) * vz[-1]) / dzi[-1]
    C[-1] += ((vz[-1] > 0) * vz[-1]) / dzi[-1]
    # vertical adection

    # shape of ni-long 1D array
    Ai[0] = -1.0 / (dzi[0]) * (Dzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / 2.0 / ysum[
        0
    ] + 1.0 / (dzi[0]) * Dzz[0] / 2.0 * (
        -1.0 / Hpi[0]
        + ms * g[0] / (Navo * kb * Ti[0])
        + alpha / Ti[0] * (Tco[1] - Tco[0]) / dzi[0]
    )
    Bi[0] = 1.0 / (dzi[0]) * (Dzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / 2.0 / ysum[
        1
    ] + 1.0 / (dzi[0]) * Dzz[0] / 2.0 * (
        -1.0 / Hpi[0]
        + ms * g[0] / (Navo * kb * Ti[0])
        + alpha / Ti[0] * (Tco[1] - Tco[0]) / dzi[0]
    )
    Ci[0] = 0
    Ai[nz - 1] = -1.0 / (dzi[-1]) * (Dzz[nz - 2] / dzi[-1]) * (
        ysum[nz - 1] + ysum[nz - 2]
    ) / 2.0 / ysum[nz - 1] - 1.0 / (dzi[-1]) * Dzz[-1] / 2.0 * (
        -1.0 / Hpi[-1]
        + ms * g[-1] / (Navo * kb * Ti[-1])
        + alpha / Ti[-1] * (Tco[-1] - Tco[-2]) / dzi[-1]
    )
    Bi[nz - 1] = 0
    Ci[nz - 1] = 1.0 / (dzi[-1]) * (Dzz[nz - 2] / dzi[-1]) * (
        ysum[nz - 1] + ysum[nz - 2]
    ) / 2.0 / ysum[nz - 2] - 1.0 / (dzi[-1]) * Dzz[-1] / 2.0 * (
        -1.0 / Hpi[-1]
        + ms * g[-1] / (Navo * kb * Ti[-1])
        + alpha / Ti[-1] * (Tco[-1] - Tco[-2]) / dzi[-1]
    )

    for j in range(1, nz - 1):
        dz_ave = 0.5 * (dzi[j - 1] + dzi[j])
        A[j] = (
            -1.0
            / dz_ave
            * (
                Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Kzz[j - 1] / dzi[j - 1] * (ysum[j] + ysum[j - 1]) / 2.0
            )
            / ysum[j]
        )
        B[j] = 1.0 / dz_ave * Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0 / ysum[j + 1]
        C[j] = (
            1.0
            / dz_ave
            * Kzz[j - 1]
            / dzi[j - 1]
            * (ysum[j] + ysum[j - 1])
            / 2.0
            / ysum[j - 1]
        )

        # vertical adection
        A[j] += -((vz[j] > 0) * vz[j] - (vz[j - 1] < 0) * vz[j - 1]) / dz_ave
        B[j] += -((vz[j] < 0) * vz[j]) / dz_ave
        C[j] += ((vz[j - 1] > 0) * vz[j - 1]) / dz_ave
        # vertical adection

        # Ai in the shape of nz*ni and Ai[j] in the shape of ni
        Ai[j] = (
            -1.0
            / dz_ave
            * (
                Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Dzz[j - 1] / dzi[j - 1] * (ysum[j] + ysum[j - 1]) / 2.0
            )
            / ysum[j]
        )
        Bi[j] = 1.0 / dz_ave * Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0 / ysum[j + 1]
        Ci[j] = (
            1.0
            / dz_ave
            * Dzz[j - 1]
            / dzi[j - 1]
            * (ysum[j] + ysum[j - 1])
            / 2.0
            / ysum[j - 1]
        )

        Ai[j] += (
            1.0
            / (2.0 * dz_ave)
            * (
                Dzz[j]
                * (
                    -1.0 / Hpi[j]
                    + ms * g[j] / (Navo * kb * Ti[j])
                    + alpha / Ti[j] * (Tco[j + 1] - Tco[j]) / dzi[j]
                )
                - Dzz[j - 1]
                * (
                    -1.0 / Hpi[j - 1]
                    + ms * g[j] / (Navo * kb * Ti[j - 1])
                    + alpha / Ti[j - 1] * (Tco[j] - Tco[j - 1]) / dzi[j - 1]
                )
            )
        )  # /ysum[j]
        Bi[j] += (
            1.0
            / (2.0 * dz_ave)
            * Dzz[j]
            * (
                -1.0 / Hpi[j]
                + ms * g[j + 1] / (Navo * kb * Ti[j])
                + alpha / Ti[j] * (Tco[j + 1] - Tco[j]) / dzi[j]
            )
        )
        Ci[j] += (
            -1.0
            / (2.0 * dz_ave)
            * Dzz[j - 1]
            * (
                -1.0 / Hpi[j - 1]
                + ms * g[j - 1] / (Navo * kb * Ti[j - 1])
                + alpha / Ti[j - 1] * (Tco[j] - Tco[j - 1]) / dzi[j - 1]
            )
        )

    tmp0 = (A[0] + Ai[0]) * y[0] + (B[0] + Bi[0]) * y[1]  # shape of ni-long 1D array
    tmp1 = np.ndarray.flatten(
        (
            np.vstack(A[1 : nz - 1]) * y[1 : (nz - 1)]
            + np.vstack(B[1 : nz - 1]) * y[1 + 1 : (nz - 1) + 1]
            + np.vstack(C[1 : nz - 1]) * y[1 - 1 : (nz - 1) - 1]
        )
    )
    tmp1 += np.ndarray.flatten(
        Ai[1 : nz - 1] * y[1 : (nz - 1)]
        + Bi[1 : nz - 1] * y[1 + 1 : (nz - 1) + 1]
        + Ci[1 : nz - 1] * y[1 - 1 : (nz - 1) - 1]
    )  # shape of (nz-2,ni)
    tmp2 = (A[nz - 1] + Ai[nz - 1]) * y[nz - 1] + (C[nz - 1] + Ci[nz - 1]) * y[nz - 2]
    diff = np.append(np.append(tmp0, tmp1), tmp2)
    diff = diff.reshape(nz, ni)

    if self.cfg.use_topflux:
        # Don't forget dz!!! -d phi/ dz
        ### the const flux has no contribution to the jacobian ###
        diff[-1] += atm.top_flux / dzi[-1]
    if self.cfg.use_botflux:
        ### the deposition term needs to be included in the jacobian!!!
        diff[0] += (atm.bot_flux - y[0] * atm.bot_vdep) / dzi[0]

    return diff

`diffdf_no_mol(y, atm)`

function of eddy diffusion without molecular diffusion, with zero-flux boundary conditions and non-uniform grids (dzi) in the form of Ajy_j + Bj+1y_j+1 + Cj-1*y_j-1

Source code in src/vulcan/op.py

def diffdf_no_mol(self, y, atm):
    """
    function of eddy diffusion without molecular diffusion, with zero-flux boundary conditions and non-uniform grids (dzi)
    in the form of Aj*y_j + Bj+1*y_j+1 + Cj-1*y_j-1
    """

    nz = self.cfg.nz

    y = y.copy()
    # TEST excluding non-gaseous species
    if self.cfg.non_gas_sp:
        ysum = np.sum(y[:, atm.gas_indx], axis=1)
    else:
        ysum = np.sum(y, axis=1)
    # TEST excluding non-gaseous species
    dzi = atm.dzi.copy()
    Kzz = atm.Kzz.copy()
    vz = atm.vz.copy()

    A, B, C = np.zeros(nz), np.zeros(nz), np.zeros(nz)

    A[0] = -1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / 2.0 / ysum[0]
    B[0] = 1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / 2.0 / ysum[1]
    C[0] = 0
    A[nz - 1] = (
        -1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[nz - 1] + ysum[nz - 2])
        / 2.0
        / ysum[nz - 1]
    )
    B[nz - 1] = 0
    C[nz - 1] = (
        1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[nz - 1] + ysum[nz - 2])
        / 2.0
        / ysum[nz - 2]
    )

    # vertical adection with zero-flux B.C.
    A[0] += -((vz[0] > 0) * vz[0]) / dzi[0]
    B[0] += -((vz[0] < 0) * vz[0]) / dzi[0]
    A[-1] += ((vz[-1] < 0) * vz[-1]) / dzi[-1]
    C[-1] += ((vz[-1] > 0) * vz[-1]) / dzi[-1]
    # vertical adection

    for j in range(1, nz - 1):
        dz_ave = 0.5 * (dzi[j - 1] + dzi[j])
        A[j] = (
            -2.0
            / (dzi[j - 1] + dzi[j])
            * (
                Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Kzz[j - 1] / dzi[j - 1] * (ysum[j] + ysum[j - 1]) / 2.0
            )
            / ysum[j]
        )
        B[j] = (
            2.0
            / (dzi[j - 1] + dzi[j])
            * Kzz[j]
            / dzi[j]
            * (ysum[j + 1] + ysum[j])
            / 2.0
            / ysum[j + 1]
        )
        C[j] = (
            2.0
            / (dzi[j - 1] + dzi[j])
            * Kzz[j - 1]
            / dzi[j - 1]
            * (ysum[j] + ysum[j - 1])
            / 2.0
            / ysum[j - 1]
        )

        # vertical adection
        A[j] += -((vz[j] > 0) * vz[j] - (vz[j - 1] < 0) * vz[j - 1]) / dz_ave
        B[j] += -((vz[j] < 0) * vz[j]) / dz_ave
        C[j] += ((vz[j - 1] > 0) * vz[j - 1]) / dz_ave
        # vertical adection

    tmp0 = A[0] * y[0] + B[0] * y[1]
    tmp1 = np.ndarray.flatten(
        (
            np.vstack(A[1 : nz - 1]) * y[1 : (nz - 1)]
            + np.vstack(B[1 : nz - 1]) * y[1 + 1 : (nz - 1) + 1]
            + np.vstack(C[1 : nz - 1]) * y[1 - 1 : (nz - 1) - 1]
        )
    )
    tmp2 = A[nz - 1] * y[nz - 1] + C[nz - 1] * y[nz - 2]
    diff = np.append(np.append(tmp0, tmp1), tmp2)
    diff = diff.reshape(nz, ni)

    if self.cfg.use_topflux:
        # Don't forget dz!!! -d phi/ dz
        ### the const flux has no contribution to the jacobian ###
        diff[-1] += atm.top_flux / dzi[-1]
    if self.cfg.use_botflux:
        ### the deposition term needs to be included in the jacobian!!!
        diff[0] += (atm.bot_flux - y[0] * atm.bot_vdep) / dzi[0]
    return diff

`diffdf_settling(y, atm)`

function of eddy diffusion including molecular diffusion and the settling velocity for particles, with zero-flux boundary conditions and non-uniform grids (dzi) in the form of Ajy_j + Bj+1y_j+1 + Cj-1*y_j-1

Source code in src/vulcan/op.py

def diffdf_settling(self, y, atm):
    """
    function of eddy diffusion including molecular diffusion and the settling velocity for particles, with zero-flux boundary conditions and non-uniform grids (dzi)
    in the form of Aj*y_j + Bj+1*y_j+1 + Cj-1*y_j-1
    """

    nz = self.cfg.nz

    y = y.copy()

    # TEST condensation excluding non-gaseous species
    if self.cfg.non_gas_sp:
        ysum = np.sum(y[:, atm.gas_indx], axis=1)
    else:
        ysum = np.sum(y, axis=1)
    # TEST condensation excluding non-gaseous species

    dzi = atm.dzi.copy()
    Kzz = atm.Kzz.copy()
    vz = atm.vz.copy()
    Dzz = atm.Dzz.copy()
    vs = atm.vs.copy()
    alpha = atm.alpha.copy()
    Tco = atm.Tco
    ms = atm.ms.copy()
    Hp = atm.Hp.copy()
    g = atm.g
    Ti = atm.Ti
    Hpi = atm.Hpi
    # # define T_1/2 for the molecular diffusion
    #         Ti = 0.5*(Tco + np.roll(Tco,-1))
    #         Ti = Ti[:-1]
    #         Hpi = 0.5*(Hp + np.roll(Hp,-1))
    #         Hpi = Hpi[:-1]
    #         # store Ti and Hpi
    #         atm.Ti = Ti
    #         atm.Hpi = Hpi

    A, B, C = np.zeros(nz), np.zeros(nz), np.zeros(nz)
    Ai, Bi, Ci = [np.zeros((nz, ni)) for i in range(3)]

    A[0] = -1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / 2.0 / ysum[0]
    B[0] = 1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / 2.0 / ysum[1]
    C[0] = 0
    A[nz - 1] = (
        -1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[nz - 1] + ysum[nz - 2])
        / 2.0
        / ysum[nz - 1]
    )
    B[nz - 1] = 0
    C[nz - 1] = (
        1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[nz - 1] + ysum[nz - 2])
        / 2.0
        / ysum[nz - 2]
    )

    # vertical adection (with closed B.C.)
    A[0] += -((vz[0] > 0) * vz[0]) / dzi[0]
    B[0] += -((vz[0] < 0) * vz[0]) / dzi[0]
    A[-1] += ((vz[-1] < 0) * vz[-1]) / dzi[-1]
    C[-1] += ((vz[-1] > 0) * vz[-1]) / dzi[-1]
    # vertical adection

    # shape of ni-long 1D array
    # Including the settling velocity of the particles
    Ai[0] = (
        -1.0 / (dzi[0]) * (Dzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / 2.0 / ysum[0]
        + 1.0
        / (dzi[0])
        * Dzz[0]
        / 2.0
        * (
            -1.0 / Hpi[0]
            + ms * g[0] / (Navo * kb * Ti[0])
            + alpha / Ti[0] * (Tco[1] - Tco[0]) / dzi[0]
        )
        - ((vs[0] > 0) * vs[0]) / dzi[0]
    )
    Bi[0] = (
        1.0 / (dzi[0]) * (Dzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / 2.0 / ysum[1]
        + 1.0
        / (dzi[0])
        * Dzz[0]
        / 2.0
        * (
            -1.0 / Hpi[0]
            + ms * g[0] / (Navo * kb * Ti[0])
            + alpha / Ti[0] * (Tco[1] - Tco[0]) / dzi[0]
        )
        - ((vs[0] < 0) * vs[0]) / dzi[0]
    )
    # Ci[0] = 0
    Ai[nz - 1] = (
        -1.0
        / (dzi[-1])
        * (Dzz[nz - 2] / dzi[-1])
        * (ysum[nz - 1] + ysum[nz - 2])
        / 2.0
        / ysum[nz - 1]
        - 1.0
        / (dzi[-1])
        * Dzz[-1]
        / 2.0
        * (
            -1.0 / Hpi[-1]
            + ms * g[-1] / (Navo * kb * Ti[-1])
            + alpha / Ti[-1] * (Tco[-1] - Tco[-2]) / dzi[-1]
        )
        + ((vs[-1] < 0) * vs[-1]) / dzi[-1]
    )
    # Bi[nz-1] = 0
    Ci[nz - 1] = (
        1.0
        / (dzi[-1])
        * (Dzz[nz - 2] / dzi[-1])
        * (ysum[nz - 1] + ysum[nz - 2])
        / 2.0
        / ysum[nz - 2]
        - 1.0
        / (dzi[-1])
        * Dzz[-1]
        / 2.0
        * (
            -1.0 / Hpi[-1]
            + ms * g[-1] / (Navo * kb * Ti[-1])
            + alpha / Ti[-1] * (Tco[-1] - Tco[-2]) / dzi[-1]
        )
        + ((vs[-1] > 0) * vs[-1]) / dzi[-1]
    )

    for j in range(1, nz - 1):
        dz_ave = 0.5 * (dzi[j - 1] + dzi[j])
        A[j] = (
            -1.0
            / dz_ave
            * (
                Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Kzz[j - 1] / dzi[j - 1] * (ysum[j] + ysum[j - 1]) / 2.0
            )
            / ysum[j]
        )
        B[j] = 1.0 / dz_ave * Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0 / ysum[j + 1]
        C[j] = (
            1.0
            / dz_ave
            * Kzz[j - 1]
            / dzi[j - 1]
            * (ysum[j] + ysum[j - 1])
            / 2.0
            / ysum[j - 1]
        )

        # vertical adection
        A[j] += -((vz[j] > 0) * vz[j] - (vz[j - 1] < 0) * vz[j - 1]) / dz_ave
        B[j] += -((vz[j] < 0) * vz[j]) / dz_ave
        C[j] += ((vz[j - 1] > 0) * vz[j - 1]) / dz_ave
        # vertical adection

        # Ai in the shape of nz*ni and Ai[j] in the shape of ni
        # Including the settling velocity of the particles
        Ai[j] = (
            -1.0
            / dz_ave
            * (
                Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Dzz[j - 1] / dzi[j - 1] * (ysum[j] + ysum[j - 1]) / 2.0
            )
            / ysum[j]
            - ((vs[j] > 0) * vs[j] - (vs[j - 1] < 0) * vs[j - 1]) / dz_ave
        )
        Bi[j] = (
            1.0 / dz_ave * Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0 / ysum[j + 1]
            - ((vs[j] < 0) * vs[j]) / dz_ave
        )
        Ci[j] = (
            1.0
            / dz_ave
            * Dzz[j - 1]
            / dzi[j - 1]
            * (ysum[j] + ysum[j - 1])
            / 2.0
            / ysum[j - 1]
            + ((vs[j - 1] > 0) * vs[j - 1]) / dz_ave
        )

        Ai[j] += (
            1.0
            / (2.0 * dz_ave)
            * (
                Dzz[j]
                * (
                    -1.0 / Hpi[j]
                    + ms * g[j] / (Navo * kb * Ti[j])
                    + alpha / Ti[j] * (Tco[j + 1] - Tco[j]) / dzi[j]
                )
                - Dzz[j - 1]
                * (
                    -1.0 / Hpi[j - 1]
                    + ms * g[j] / (Navo * kb * Ti[j - 1])
                    + alpha / Ti[j - 1] * (Tco[j] - Tco[j - 1]) / dzi[j - 1]
                )
            )
        )  # /ysum[j]
        Bi[j] += (
            1.0
            / (2.0 * dz_ave)
            * Dzz[j]
            * (
                -1.0 / Hpi[j]
                + ms * g[j + 1] / (Navo * kb * Ti[j])
                + alpha / Ti[j] * (Tco[j + 1] - Tco[j]) / dzi[j]
            )
        )
        Ci[j] += (
            -1.0
            / (2.0 * dz_ave)
            * Dzz[j - 1]
            * (
                -1.0 / Hpi[j - 1]
                + ms * g[j - 1] / (Navo * kb * Ti[j - 1])
                + alpha / Ti[j - 1] * (Tco[j] - Tco[j - 1]) / dzi[j - 1]
            )
        )

    tmp0 = (A[0] + Ai[0]) * y[0] + (B[0] + Bi[0]) * y[1]  # shape of ni-long 1D array
    tmp1 = np.ndarray.flatten(
        (
            np.vstack(A[1 : nz - 1]) * y[1 : (nz - 1)]
            + np.vstack(B[1 : nz - 1]) * y[1 + 1 : (nz - 1) + 1]
            + np.vstack(C[1 : nz - 1]) * y[1 - 1 : (nz - 1) - 1]
        )
    )
    tmp1 += np.ndarray.flatten(
        Ai[1 : nz - 1] * y[1 : (nz - 1)]
        + Bi[1 : nz - 1] * y[1 + 1 : (nz - 1) + 1]
        + Ci[1 : nz - 1] * y[1 - 1 : (nz - 1) - 1]
    )  # shape of (nz-2,ni)
    tmp2 = (A[nz - 1] + Ai[nz - 1]) * y[nz - 1] + (C[nz - 1] + Ci[nz - 1]) * y[nz - 2]
    diff = np.append(np.append(tmp0, tmp1), tmp2)
    diff = diff.reshape(nz, ni)

    if self.cfg.use_topflux:
        # Don't forget dz!!! -d phi/ dz
        ### the const flux has no contribution to the jacobian ###
        diff[-1] += atm.top_flux / dzi[-1]
    if self.cfg.use_botflux:
        ### the deposition term needs to be included in the jacobian!!!
        diff[0] += (atm.bot_flux - y[0] * atm.bot_vdep) / dzi[0]

    return diff

`diffdf_settling_vm(y, atm)`

added vm for molecular diffusion function of eddy diffusion including molecular diffusion and the settling velocity for particles, with zero-flux boundary conditions and non-uniform grids (dzi) in the form of Ajy_j + Bj+1y_j+1 + Cj-1*y_j-1

Source code in src/vulcan/op.py

def diffdf_settling_vm(self, y, atm):
    """
    added vm for molecular diffusion
    function of eddy diffusion including molecular diffusion and the settling velocity for particles, with zero-flux boundary conditions and non-uniform grids (dzi)
    in the form of Aj*y_j + Bj+1*y_j+1 + Cj-1*y_j-1
    """

    nz = self.cfg.nz
    y = y.copy()

    if self.cfg.non_gas_sp:
        ysum = np.sum(y[:, atm.gas_indx], axis=1)
    else:
        ysum = np.sum(y, axis=1)

    dzi = atm.dzi.copy()
    Kzz = atm.Kzz.copy()
    vz = atm.vz.copy()
    Dzz = atm.Dzz.copy()
    vs = atm.vs.copy()
    alpha = atm.alpha.copy()
    Tco = atm.Tco.copy()
    ms = atm.ms.copy()
    Hp = atm.Hp.copy()
    g = atm.g
    Ti = atm.Ti
    Hpi = atm.Hpi

    vm = atm.vm
    # shape: nz x ni
    # vm defined in build.py
    # vm = - Dzz_cen * ( ms[np.newaxis,:]*g[:,np.newaxis]/(Navo*kb*Tco[:,np.newaxis]) - 1./Hp[:,np.newaxis] +  alpha/Tco[:,np.newaxis]*(delta_T[:,np.newaxis])/dz[:,np.newaxis]  )

    A, B, C = np.zeros(nz), np.zeros(nz), np.zeros(nz)
    Ai, Bi, Ci = [np.zeros((nz, ni)) for i in range(3)]

    A[0] = -1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / 2.0 / ysum[0]
    B[0] = 1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / 2.0 / ysum[1]
    C[0] = 0
    A[nz - 1] = (
        -1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[nz - 1] + ysum[nz - 2])
        / 2.0
        / ysum[nz - 1]
    )
    B[nz - 1] = 0
    C[nz - 1] = (
        1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[nz - 1] + ysum[nz - 2])
        / 2.0
        / ysum[nz - 2]
    )

    # vertical adection (with closed B.C.)
    A[0] += -((vz[0] > 0) * vz[0]) / dzi[0]
    B[0] += -((vz[0] < 0) * vz[0]) / dzi[0]
    A[-1] += ((vz[-1] < 0) * vz[-1]) / dzi[-1]
    C[-1] += ((vz[-1] > 0) * vz[-1]) / dzi[-1]
    # vertical adection

    # shape of ni-long 1D array
    # Including the settling velocity of the particles and the advective component of molecular diffusion
    Ai[0] = (
        -1.0 / (dzi[0]) * (Dzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / 2.0 / ysum[0]
        - ((vs[0] > 0) * vs[0]) / dzi[0]
    )
    Bi[0] = (
        1.0 / (dzi[0]) * (Dzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / 2.0 / ysum[1]
        - ((vs[0] < 0) * vs[0]) / dzi[0]
    )
    # Ci[0] = 0
    Ai[nz - 1] = (
        -1.0
        / (dzi[-1])
        * (Dzz[nz - 2] / dzi[-1])
        * (ysum[nz - 1] + ysum[nz - 2])
        / 2.0
        / ysum[nz - 1]
        + ((vm[-1] < 0) * vm[-1]) / dzi[-1]
        + ((vs[-1] < 0) * vs[-1]) / dzi[-1]
    )
    # Bi[nz-1] = 0
    Ci[nz - 1] = (
        1.0
        / (dzi[-1])
        * (Dzz[nz - 2] / dzi[-1])
        * (ysum[nz - 1] + ysum[nz - 2])
        / 2.0
        / ysum[nz - 2]
        + ((vm[-1] > 0) * vm[-1]) / dzi[-1]
        + ((vs[-1] > 0) * vs[-1]) / dzi[-1]
    )

    for j in range(1, nz - 1):
        dz_ave = 0.5 * (dzi[j - 1] + dzi[j])
        A[j] = (
            -1.0
            / dz_ave
            * (
                Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Kzz[j - 1] / dzi[j - 1] * (ysum[j] + ysum[j - 1]) / 2.0
            )
            / ysum[j]
        )
        B[j] = 1.0 / dz_ave * Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0 / ysum[j + 1]
        C[j] = (
            1.0
            / dz_ave
            * Kzz[j - 1]
            / dzi[j - 1]
            * (ysum[j] + ysum[j - 1])
            / 2.0
            / ysum[j - 1]
        )

        # vertical adection
        A[j] += -((vz[j] > 0) * vz[j] - (vz[j - 1] < 0) * vz[j - 1]) / dz_ave
        B[j] += -((vz[j] < 0) * vz[j]) / dz_ave
        C[j] += ((vz[j - 1] > 0) * vz[j - 1]) / dz_ave
        # vertical adection

        # Ai in the shape of nz*ni and Ai[j] in the shape of ni
        # Including the settling velocity of the particles

        # diffusion component
        Ai[j] = (
            -1.0
            / dz_ave
            * (
                Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Dzz[j - 1] / dzi[j - 1] * (ysum[j] + ysum[j - 1]) / 2.0
            )
            / ysum[j]
        )
        Bi[j] = 1.0 / dz_ave * Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0 / ysum[j + 1]
        Ci[j] = (
            1.0
            / dz_ave
            * Dzz[j - 1]
            / dzi[j - 1]
            * (ysum[j] + ysum[j - 1])
            / 2.0
            / ysum[j - 1]
        )
        # diffusion component

        # advective component using upwind (inc. from Dzz and from vs)
        Ai[j] += (
            -((vm[j] > 0) * vm[j] - (vm[j - 1] < 0) * vm[j - 1]) / dz_ave
            - ((vs[j] > 0) * vs[j] - (vs[j - 1] < 0) * vs[j - 1]) / dz_ave
        )
        Bi[j] += -((vm[j] < 0) * vm[j]) / dz_ave - ((vs[j] < 0) * vs[j]) / dz_ave
        Ci[j] += (
            +((vm[j - 1] > 0) * vm[j - 1]) / dz_ave + ((vs[j - 1] > 0) * vs[j - 1]) / dz_ave
        )
        # advective component using upwind

    tmp0 = (A[0] + Ai[0]) * y[0] + (B[0] + Bi[0]) * y[1]  # shape of ni-long 1D array
    tmp1 = np.ndarray.flatten(
        (
            np.vstack(A[1 : nz - 1]) * y[1 : (nz - 1)]
            + np.vstack(B[1 : nz - 1]) * y[1 + 1 : (nz - 1) + 1]
            + np.vstack(C[1 : nz - 1]) * y[1 - 1 : (nz - 1) - 1]
        )
    )
    tmp1 += np.ndarray.flatten(
        Ai[1 : nz - 1] * y[1 : (nz - 1)]
        + Bi[1 : nz - 1] * y[1 + 1 : (nz - 1) + 1]
        + Ci[1 : nz - 1] * y[1 - 1 : (nz - 1) - 1]
    )  # shape of (nz-2,ni)
    tmp2 = (A[nz - 1] + Ai[nz - 1]) * y[nz - 1] + (C[nz - 1] + Ci[nz - 1]) * y[nz - 2]
    diff = np.append(np.append(tmp0, tmp1), tmp2)
    diff = diff.reshape(nz, ni)

    if self.cfg.use_topflux:
        # Don't forget dz!!! -d phi/ dz
        ### the const flux has no contribution to the jacobian ###
        diff[-1] += atm.top_flux / dzi[-1]
    if self.cfg.use_botflux:
        ### the deposition term needs to be included in the jacobian!!!
        diff[0] += (atm.bot_flux - y[0] * atm.bot_vdep) / dzi[0]

    return diff

`jac_tot(var, atm)`

jacobian matrix for dn/dt + dphi/dz = P - L (including molecular diffusion) zero-flux BC: 1st derivitive of y is zero

Source code in src/vulcan/op.py

def jac_tot(self, var, atm):
    """
    jacobian matrix for dn/dt + dphi/dz = P - L (including molecular diffusion)
    zero-flux BC:  1st derivitive of y is zero
    """

    nz = self.cfg.nz
    y = var.y.copy()

    # TEST condensation excluding non-gaseous species
    if self.cfg.non_gas_sp:
        ysum = np.sum(y[:, atm.gas_indx], axis=1)
    else:
        ysum = np.sum(y, axis=1)
    # TEST condensation excluding non-gaseous species
    dzi = atm.dzi.copy()
    Kzz = atm.Kzz.copy()
    Dzz = atm.Dzz.copy()
    vz = atm.vz.copy()
    alpha = atm.alpha.copy()
    Tco = atm.Tco
    mu, ms = atm.mu.copy(), atm.ms.copy()
    g = atm.g

    # define T_1/2 for the molecular diffusion
    # Ti = 0.5*(Tco + np.roll(Tco,-1))
    # Ti = Ti[:-1]

    Ti = atm.Ti.copy()
    Hpi = atm.Hpi.copy()

    dfdy = achemjac(y, atm.M, var.k)
    j_indx = []

    for j in range(nz):
        j_indx.append(np.arange(j * ni, j * ni + ni))

    for j in range(1, nz - 1):
        # excluding the buttom and the top cell
        # at j level consists of ni species
        dz_ave = 0.5 * (dzi[j - 1] + dzi[j])
        dfdy[j_indx[j], j_indx[j]] += (
            -1.0
            / dz_ave
            * (
                Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Kzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / 2.0
            )
            / ysum[j]
            - ((vz[j] > 0) * vz[j] - (vz[j - 1] < 0) * vz[j - 1]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j + 1]] += (
            1.0 / dz_ave * (Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / (2.0 * ysum[j + 1]))
            - ((vz[j] < 0) * vz[j]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j - 1]] += (
            1.0
            / dz_ave
            * (Kzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / (2.0 * ysum[j - 1]))
            + ((vz[j - 1] > 0) * vz[j - 1]) / dz_ave
        )

        # [j_indx[j], j_indx[j]] has size ni*ni
        dfdy[j_indx[j], j_indx[j]] += -1.0 / dz_ave * (
            Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
            + Dzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / 2.0
        ) / ysum[j] + 1.0 / (2.0 * dz_ave) * (
            Dzz[j]
            * (
                -1.0 / Hpi[j]
                + ms * g[j] / (Navo * kb * Ti[j])
                + alpha / Ti[j] * (Tco[j + 1] - Tco[j]) / dzi[j]
            )
            - Dzz[j - 1]
            * (
                -1.0 / Hpi[j - 1]
                + ms * g[j] / (Navo * kb * Ti[j - 1])
                + alpha / Ti[j - 1] * (Tco[j] - Tco[j - 1]) / dzi[j - 1]
            )
        )
        dfdy[j_indx[j], j_indx[j + 1]] += 1.0 / dz_ave * (
            Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / (2.0 * ysum[j + 1])
        ) + 1.0 / (2.0 * dz_ave) * Dzz[j] * (
            -1.0 / Hpi[j]
            + ms * g[j + 1] / (Navo * kb * Ti[j])
            + alpha / Ti[j] * (Tco[j + 1] - Tco[j]) / dzi[j]
        )
        dfdy[j_indx[j], j_indx[j - 1]] += 1.0 / dz_ave * (
            Dzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / (2.0 * ysum[j - 1])
        ) - 1.0 / (2.0 * dz_ave) * Dzz[j - 1] * (
            -1.0 / Hpi[j - 1]
            + ms * g[j - 1] / (Navo * kb * Ti[j - 1])
            + alpha / Ti[j - 1] * (Tco[j] - Tco[j - 1]) / dzi[j - 1]
        )

    dfdy[j_indx[0], j_indx[0]] += (
        -1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[0])
        - ((vz[0] > 0) * vz[0]) / dzi[0]
    )
    dfdy[j_indx[0], j_indx[0]] += -1.0 / (dzi[0]) * (Dzz[0] / dzi[0]) * (
        ysum[1] + ysum[0]
    ) / (2.0 * ysum[0]) + 1.0 / (dzi[0]) * Dzz[0] / 2.0 * (
        -1.0 / Hpi[0]
        + ms * g[0] / (Navo * kb * Ti[0])
        + alpha / Ti[0] * (Tco[1] - Tco[0]) / dzi[0]
    )
    # deposition velocity
    if self.cfg.use_botflux:
        dfdy[j_indx[0], j_indx[0]] += -1.0 * atm.bot_vdep / dzi[0]

    dfdy[j_indx[0], j_indx[1]] += (
        1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[1])
        - ((vz[0] < 0) * vz[0]) / dzi[0]
    )
    dfdy[j_indx[0], j_indx[1]] += 1.0 / (dzi[0]) * (Dzz[0] / dzi[0]) * (
        ysum[1] + ysum[0]
    ) / (2.0 * ysum[1]) + 1.0 / (dzi[0]) * Dzz[0] / 2.0 * (
        -1.0 / Hpi[0]
        + ms * g[0] / (Navo * kb * Ti[0])
        + alpha / Ti[0] * (Tco[1] - Tco[0]) / dzi[0]
    )

    dfdy[j_indx[nz - 1], j_indx[nz - 1]] += (
        -1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[(nz - 1) - 1] + ysum[nz - 1])
        / (2.0 * ysum[nz - 1])
        + ((vz[-1] < 0) * vz[-1]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[nz - 1]] += -1.0 / (dzi[nz - 2]) * (
        Dzz[nz - 2] / dzi[nz - 2]
    ) * (ysum[nz - 1] + ysum[nz - 2]) / (2.0 * ysum[nz - 1]) - 1.0 / (dzi[-1]) * Dzz[
        -1
    ] / 2.0 * (
        -1.0 / Hpi[-1]
        + ms * g[-1] / (Navo * kb * Ti[-1])
        + alpha / Ti[-1] * (Tco[-1] - Tco[-2]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[(nz - 1) - 1]] += (
        1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[(nz - 1) - 1] + ysum[nz - 1])
        / (2.0 * ysum[(nz - 1) - 1])
        + ((vz[-1] > 0) * vz[-1]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[(nz - 1) - 1]] += 1.0 / (dzi[nz - 2]) * (
        Dzz[nz - 2] / dzi[nz - 2]
    ) * (ysum[nz - 1] + ysum[nz - 2]) / (2.0 * ysum[(nz - 1) - 1]) - 1.0 / (dzi[-1]) * Dzz[
        -1
    ] / 2.0 * (
        -1.0 / Hpi[-1]
        + ms * g[-1] / (Navo * kb * Ti[-1])
        + alpha / Ti[-1] * (Tco[-1] - Tco[-2]) / dzi[-1]
    )

    return dfdy

`lhs_jac_fix_all_bot(var, atm)`

directly constructing lhs = 1./(rh)sparse.identity(ni*nz) - dfdy jacobian matrix for dn/dt + dphi/dz = P - L (including molecular diffusion) Fixed all species BC: all species at bottom (y[0]) remains fixed

Source code in src/vulcan/op.py

def lhs_jac_fix_all_bot(self, var, atm):
    """
    directly constructing lhs = 1./(r*h)*sparse.identity(ni*nz) - dfdy
    jacobian matrix for dn/dt + dphi/dz = P - L (including molecular diffusion)
    Fixed all species BC: all species at bottom (y[0]) remains fixed
    """

    nz = self.cfg.nz

    y = var.y.copy()
    # TEST condensation excluding non-gaseous species
    if self.cfg.non_gas_sp:
        ysum = np.sum(y[:, atm.gas_indx], axis=1)
    else:
        ysum = np.sum(y, axis=1)
    # TEST condensation excluding non-gaseous species
    dzi = atm.dzi.copy()
    Kzz = atm.Kzz.copy()
    Dzz = atm.Dzz.copy()
    vz = atm.vz.copy()
    alpha = atm.alpha.copy()
    Tco = atm.Tco
    mu, ms = atm.mu.copy(), atm.ms.copy()
    g = atm.g

    Ti = atm.Ti.copy()
    Hpi = atm.Hpi.copy()

    r = 1.0 + 1.0 / 2.0**0.5
    c0 = 1.0 / (r * var.dt)
    dfdy = neg_achemjac(y, atm.M, var.k, self.cfg.nz)
    np.fill_diagonal(dfdy, c0 + np.diag(dfdy))
    j_indx = []

    for j in range(nz):
        j_indx.append(np.arange(j * ni, j * ni + ni))

    for j in range(1, nz - 1):
        # excluding the buttom and the top cell
        # at j level consists of ni species
        dz_ave = 0.5 * (dzi[j - 1] + dzi[j])
        dfdy[j_indx[j], j_indx[j]] -= (
            -1.0
            / dz_ave
            * (
                Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Kzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / 2.0
            )
            / ysum[j]
            - ((vz[j] > 0) * vz[j] - (vz[j - 1] < 0) * vz[j - 1]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j + 1]] -= (
            1.0 / dz_ave * (Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / (2.0 * ysum[j + 1]))
            - ((vz[j] < 0) * vz[j]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j - 1]] -= (
            1.0
            / dz_ave
            * (Kzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / (2.0 * ysum[j - 1]))
            + ((vz[j - 1] > 0) * vz[j - 1]) / dz_ave
        )

        # [j_indx[j], j_indx[j]] has size ni*ni
        dfdy[j_indx[j], j_indx[j]] -= -1.0 / dz_ave * (
            Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
            + Dzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / 2.0
        ) / ysum[j] + 1.0 / (2.0 * dz_ave) * (
            Dzz[j]
            * (
                -1.0 / Hpi[j]
                + ms * g[j] / (Navo * kb * Ti[j])
                + alpha / Ti[j] * (Tco[j + 1] - Tco[j]) / dzi[j]
            )
            - Dzz[j - 1]
            * (
                -1.0 / Hpi[j - 1]
                + ms * g[j] / (Navo * kb * Ti[j - 1])
                + alpha / Ti[j - 1] * (Tco[j] - Tco[j - 1]) / dzi[j - 1]
            )
        )
        dfdy[j_indx[j], j_indx[j + 1]] -= 1.0 / dz_ave * (
            Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / (2.0 * ysum[j + 1])
        ) + 1.0 / (2.0 * dz_ave) * Dzz[j] * (
            -1.0 / Hpi[j]
            + ms * g[j + 1] / (Navo * kb * Ti[j])
            + alpha / Ti[j] * (Tco[j + 1] - Tco[j]) / dzi[j]
        )
        dfdy[j_indx[j], j_indx[j - 1]] -= 1.0 / dz_ave * (
            Dzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / (2.0 * ysum[j - 1])
        ) - 1.0 / (2.0 * dz_ave) * Dzz[j - 1] * (
            -1.0 / Hpi[j - 1]
            + ms * g[j - 1] / (Navo * kb * Ti[j - 1])
            + alpha / Ti[j - 1] * (Tco[j] - Tco[j - 1]) / dzi[j - 1]
        )

    # deposition velocity (off with fixed all BC)
    # if self.cfg.use_botflux : dfdy[j_indx[0], j_indx[0]] -= -1.*atm.bot_vdep /dzi[0]

    # diffusion-limited escape
    if self.cfg.diff_esc:  # not empty list
        diff_lim = np.zeros(ni)
        for sp in self.cfg.diff_esc:
            if y[-1, species.index(sp)] > 0:
                diff_lim[species.index(sp)] += (
                    atm.top_flux[species.index(sp)] / y[-1, species.index(sp)]
                )
        dfdy[j_indx[-1], j_indx[-1]] -= diff_lim  # negative

    # Fix bottom BC
    dfdy[:, j_indx[0]] = 0.0

    dfdy[j_indx[0], j_indx[1]] -= (
        1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[1])
        - ((vz[0] < 0) * vz[0]) / dzi[0]
    )
    dfdy[j_indx[0], j_indx[1]] -= 1.0 / (dzi[0]) * (Dzz[0] / dzi[0]) * (
        ysum[1] + ysum[0]
    ) / (2.0 * ysum[1]) + 1.0 / (dzi[0]) * Dzz[0] / 2.0 * (
        -1.0 / Hpi[0]
        + ms * g[0] / (Navo * kb * Ti[0])
        + alpha / Ti[0] * (Tco[1] - Tco[0]) / dzi[0]
    )

    dfdy[j_indx[nz - 1], j_indx[nz - 1]] -= (
        -1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[(nz - 1) - 1] + ysum[nz - 1])
        / (2.0 * ysum[nz - 1])
        + ((vz[-1] < 0) * vz[-1]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[nz - 1]] -= -1.0 / (dzi[nz - 2]) * (
        Dzz[nz - 2] / dzi[nz - 2]
    ) * (ysum[nz - 1] + ysum[nz - 2]) / (2.0 * ysum[nz - 1]) - 1.0 / (dzi[-1]) * Dzz[
        -1
    ] / 2.0 * (
        -1.0 / Hpi[-1]
        + ms * g[-1] / (Navo * kb * Ti[-1])
        + alpha / Ti[-1] * (Tco[-1] - Tco[-2]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[(nz - 1) - 1]] -= (
        1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[(nz - 1) - 1] + ysum[nz - 1])
        / (2.0 * ysum[(nz - 1) - 1])
        + ((vz[-1] > 0) * vz[-1]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[(nz - 1) - 1]] -= 1.0 / (dzi[nz - 2]) * (
        Dzz[nz - 2] / dzi[nz - 2]
    ) * (ysum[nz - 1] + ysum[nz - 2]) / (2.0 * ysum[(nz - 1) - 1]) - 1.0 / (dzi[-1]) * Dzz[
        -1
    ] / 2.0 * (
        -1.0 / Hpi[-1]
        + ms * g[-1] / (Navo * kb * Ti[-1])
        + alpha / Ti[-1] * (Tco[-1] - Tco[-2]) / dzi[-1]
    )

    return dfdy

`lhs_jac_no_mol(var, atm)`

directly constructing lhs = 1./(rh)sparse.identity(ni*nz) - dfdy jacobian matrix for dn/dt + dphi/dz = P - L (WITHOUT molecular diffusion) zero-flux BC: 1st derivitive of y is zero

Source code in src/vulcan/op.py

def lhs_jac_no_mol(self, var, atm):
    """
    directly constructing lhs = 1./(r*h)*sparse.identity(ni*nz) - dfdy
    jacobian matrix for dn/dt + dphi/dz = P - L (WITHOUT molecular diffusion)
    zero-flux BC:  1st derivitive of y is zero
    """

    nz = self.cfg.nz

    y = var.y.copy()
    # TEST condensation excluding non-gaseous species
    if self.cfg.non_gas_sp:
        ysum = np.sum(y[:, atm.gas_indx], axis=1)
    else:
        ysum = np.sum(y, axis=1)
    # TEST condensation excluding non-gaseous species
    dzi = atm.dzi.copy()
    Kzz = atm.Kzz.copy()
    vz = atm.vz.copy()
    Tco = atm.Tco
    mu, ms = atm.mu.copy(), atm.ms.copy()

    r = 1.0 + 1.0 / 2.0**0.5
    c0 = 1.0 / (r * var.dt)
    dfdy = neg_achemjac(y, atm.M, var.k, self.cfg.nz)
    np.fill_diagonal(dfdy, c0 + np.diag(dfdy))
    j_indx = []

    for j in range(nz):
        j_indx.append(np.arange(j * ni, j * ni + ni))

    for j in range(1, nz - 1):
        # excluding the buttom and the top cell
        # at j level consists of ni species
        dz_ave = 0.5 * (dzi[j - 1] + dzi[j])
        dfdy[j_indx[j], j_indx[j]] -= (
            -1.0
            / dz_ave
            * (
                Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Kzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / 2.0
            )
            / ysum[j]
            - ((vz[j] > 0) * vz[j] - (vz[j - 1] < 0) * vz[j - 1]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j + 1]] -= (
            1.0 / dz_ave * (Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / (2.0 * ysum[j + 1]))
            - ((vz[j] < 0) * vz[j]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j - 1]] -= (
            1.0
            / dz_ave
            * (Kzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / (2.0 * ysum[j - 1]))
            + ((vz[j - 1] > 0) * vz[j - 1]) / dz_ave
        )

    dfdy[j_indx[0], j_indx[0]] -= (
        -1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[0])
        - ((vz[0] > 0) * vz[0]) / dzi[0]
    )
    # deposition velocity
    if self.cfg.use_botflux:
        dfdy[j_indx[0], j_indx[0]] -= -1.0 * atm.bot_vdep / dzi[0]

    dfdy[j_indx[0], j_indx[1]] -= (
        1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[1])
        - ((vz[0] < 0) * vz[0]) / dzi[0]
    )

    dfdy[j_indx[nz - 1], j_indx[nz - 1]] -= (
        -1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[(nz - 1) - 1] + ysum[nz - 1])
        / (2.0 * ysum[nz - 1])
        + ((vz[-1] < 0) * vz[-1]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[(nz - 1) - 1]] -= (
        1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[(nz - 1) - 1] + ysum[nz - 1])
        / (2.0 * ysum[(nz - 1) - 1])
        + ((vz[-1] > 0) * vz[-1]) / dzi[-1]
    )

    return dfdy

`lhs_jac_no_mol_fix_all_bot(var, atm)`

directly constructing lhs = 1./(rh)sparse.identity(ni*nz) - dfdy jacobian matrix for dn/dt + dphi/dz = P - L (WITHOUT molecular diffusion) Fixed all species BC: all species at bottom (y[0]) remains fixed

Source code in src/vulcan/op.py

def lhs_jac_no_mol_fix_all_bot(self, var, atm):
    """
    directly constructing lhs = 1./(r*h)*sparse.identity(ni*nz) - dfdy
    jacobian matrix for dn/dt + dphi/dz = P - L (WITHOUT molecular diffusion)
    Fixed all species BC: all species at bottom (y[0]) remains fixed
    """

    nz = self.cfg.nz

    y = var.y.copy()
    # TEST condensation excluding non-gaseous species
    if self.cfg.non_gas_sp:
        ysum = np.sum(y[:, atm.gas_indx], axis=1)
    else:
        ysum = np.sum(y, axis=1)
    # TEST condensation excluding non-gaseous species
    dzi = atm.dzi.copy()
    Kzz = atm.Kzz.copy()
    vz = atm.vz.copy()
    Tco = atm.Tco
    mu, ms = atm.mu.copy(), atm.ms.copy()

    r = 1.0 + 1.0 / 2.0**0.5
    c0 = 1.0 / (r * var.dt)
    dfdy = neg_achemjac(y, atm.M, var.k, self.cfg.nz)
    np.fill_diagonal(dfdy, c0 + np.diag(dfdy))
    j_indx = []

    for j in range(nz):
        j_indx.append(np.arange(j * ni, j * ni + ni))

    for j in range(1, nz - 1):
        # excluding the buttom and the top cell
        # at j level consists of ni species
        dz_ave = 0.5 * (dzi[j - 1] + dzi[j])
        dfdy[j_indx[j], j_indx[j]] -= (
            -1.0
            / dz_ave
            * (
                Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Kzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / 2.0
            )
            / ysum[j]
            - ((vz[j] > 0) * vz[j] - (vz[j - 1] < 0) * vz[j - 1]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j + 1]] -= (
            1.0 / dz_ave * (Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / (2.0 * ysum[j + 1]))
            - ((vz[j] < 0) * vz[j]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j - 1]] -= (
            1.0
            / dz_ave
            * (Kzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / (2.0 * ysum[j - 1]))
            + ((vz[j - 1] > 0) * vz[j - 1]) / dz_ave
        )

    # dfdy[j_indx[0], j_indx[0]] -= -1./(dzi[0])*(Kzz[0]/dzi[0]) * (ysum[1]+ysum[0])/(2.*ysum[0]) -( (vz[0]>0)*vz[0] )/dzi[0]
    # deposition velocity (off with fixed all BC)
    # if self.cfg.use_botflux : dfdy[j_indx[0], j_indx[0]] -= -1.*atm.bot_vdep /dzi[0]

    # Fix bottom BC
    dfdy[:, j_indx[0]] = 0.0

    dfdy[j_indx[0], j_indx[1]] -= (
        1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[1])
        - ((vz[0] < 0) * vz[0]) / dzi[0]
    )

    dfdy[j_indx[nz - 1], j_indx[nz - 1]] -= (
        -1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[(nz - 1) - 1] + ysum[nz - 1])
        / (2.0 * ysum[nz - 1])
        + ((vz[-1] < 0) * vz[-1]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[(nz - 1) - 1]] -= (
        1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[(nz - 1) - 1] + ysum[nz - 1])
        / (2.0 * ysum[(nz - 1) - 1])
        + ((vz[-1] > 0) * vz[-1]) / dzi[-1]
    )

    return dfdy

`lhs_jac_settling(var, atm)`

directly constructing lhs = 1./(rh)sparse.identity(ni*nz) - dfdy jacobian matrix for dn/dt + dphi/dz = P - L (including molecular diffusion and gravitation settling for particles) zero-flux BC: 1st derivitive of y is zero

Source code in src/vulcan/op.py

def lhs_jac_settling(self, var, atm):
    """
    directly constructing lhs = 1./(r*h)*sparse.identity(ni*nz) - dfdy
    jacobian matrix for dn/dt + dphi/dz = P - L (including molecular diffusion and gravitation settling for particles)
    zero-flux BC:  1st derivitive of y is zero
    """

    nz = self.cfg.nz

    y = var.y.copy()
    # TEST condensation excluding non-gaseous species
    if self.cfg.non_gas_sp:
        ysum = np.sum(y[:, atm.gas_indx], axis=1)
    else:
        ysum = np.sum(y, axis=1)
    # TEST condensation excluding non-gaseous species
    dzi = atm.dzi.copy()
    Kzz = atm.Kzz.copy()
    Dzz = atm.Dzz.copy()
    vz = atm.vz.copy()
    vs = atm.vs.copy()
    alpha = atm.alpha.copy()
    Tco = atm.Tco
    mu, ms = atm.mu.copy(), atm.ms.copy()
    g = atm.g

    Ti = atm.Ti.copy()
    Hpi = atm.Hpi.copy()

    # c0 = 1./(r*h) where r = 1. + 1./2.**0.5
    r = 1.0 + 1.0 / 2.0**0.5
    c0 = 1.0 / (r * var.dt)
    dfdy = neg_achemjac(y, atm.M, var.k, self.cfg.nz)
    np.fill_diagonal(dfdy, c0 + np.diag(dfdy))
    j_indx = []

    for j in range(nz):
        j_indx.append(np.arange(j * ni, j * ni + ni))

    for j in range(1, nz - 1):
        # excluding the buttom and the top cell
        # at j level consists of ni species
        dz_ave = 0.5 * (dzi[j - 1] + dzi[j])
        dfdy[j_indx[j], j_indx[j]] -= (
            -1.0
            / dz_ave
            * (
                Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Kzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / 2.0
            )
            / ysum[j]
            - ((vz[j] > 0) * vz[j] - (vz[j - 1] < 0) * vz[j - 1]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j + 1]] -= (
            1.0 / dz_ave * (Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / (2.0 * ysum[j + 1]))
            - ((vz[j] < 0) * vz[j]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j - 1]] -= (
            1.0
            / dz_ave
            * (Kzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / (2.0 * ysum[j - 1]))
            + ((vz[j - 1] > 0) * vz[j - 1]) / dz_ave
        )

        # [j_indx[j], j_indx[j]] has size ni*ni
        dfdy[j_indx[j], j_indx[j]] -= (
            -1.0
            / dz_ave
            * (
                Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Dzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / 2.0
            )
            / ysum[j]
            + 1.0
            / (2.0 * dz_ave)
            * (
                Dzz[j]
                * (
                    -1.0 / Hpi[j]
                    + ms * g[j] / (Navo * kb * Ti[j])
                    + alpha / Ti[j] * (Tco[j + 1] - Tco[j]) / dzi[j]
                )
                - Dzz[j - 1]
                * (
                    -1.0 / Hpi[j - 1]
                    + ms * g[j] / (Navo * kb * Ti[j - 1])
                    + alpha / Ti[j - 1] * (Tco[j] - Tco[j - 1]) / dzi[j - 1]
                )
            )
            - ((vs[j] > 0) * vs[j] - (vs[j - 1] < 0) * vs[j - 1]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j + 1]] -= (
            1.0 / dz_ave * (Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / (2.0 * ysum[j + 1]))
            + 1.0
            / (2.0 * dz_ave)
            * Dzz[j]
            * (
                -1.0 / Hpi[j]
                + ms * g[j + 1] / (Navo * kb * Ti[j])
                + alpha / Ti[j] * (Tco[j + 1] - Tco[j]) / dzi[j]
            )
            - ((vs[j] < 0) * vs[j]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j - 1]] -= (
            1.0
            / dz_ave
            * (Dzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / (2.0 * ysum[j - 1]))
            - 1.0
            / (2.0 * dz_ave)
            * Dzz[j - 1]
            * (
                -1.0 / Hpi[j - 1]
                + ms * g[j - 1] / (Navo * kb * Ti[j - 1])
                + alpha / Ti[j - 1] * (Tco[j] - Tco[j - 1]) / dzi[j - 1]
            )
            + ((vs[j - 1] > 0) * vs[j - 1]) / dz_ave
        )

    dfdy[j_indx[0], j_indx[0]] -= (
        -1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[0])
        - ((vz[0] > 0) * vz[0]) / dzi[0]
    )
    dfdy[j_indx[0], j_indx[0]] -= (
        -1.0 / (dzi[0]) * (Dzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[0])
        + 1.0
        / (dzi[0])
        * Dzz[0]
        / 2.0
        * (
            -1.0 / Hpi[0]
            + ms * g[0] / (Navo * kb * Ti[0])
            + alpha / Ti[0] * (Tco[1] - Tco[0]) / dzi[0]
        )
        - ((vs[0] > 0) * vs[0]) / dzi[0]
    )
    # deposition velocity
    if self.cfg.use_botflux:
        dfdy[j_indx[0], j_indx[0]] -= -1.0 * atm.bot_vdep / dzi[0]

    dfdy[j_indx[0], j_indx[1]] -= (
        1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[1])
        - ((vz[0] < 0) * vz[0]) / dzi[0]
    )
    dfdy[j_indx[0], j_indx[1]] -= (
        1.0 / (dzi[0]) * (Dzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[1])
        + 1.0
        / (dzi[0])
        * Dzz[0]
        / 2.0
        * (
            -1.0 / Hpi[0]
            + ms * g[0] / (Navo * kb * Ti[0])
            + alpha / Ti[0] * (Tco[1] - Tco[0]) / dzi[0]
        )
        - ((vs[0] < 0) * vs[0]) / dzi[0]
    )

    dfdy[j_indx[nz - 1], j_indx[nz - 1]] -= (
        -1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[(nz - 1) - 1] + ysum[nz - 1])
        / (2.0 * ysum[nz - 1])
        + ((vz[-1] < 0) * vz[-1]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[nz - 1]] -= (
        -1.0
        / (dzi[nz - 2])
        * (Dzz[nz - 2] / dzi[nz - 2])
        * (ysum[nz - 1] + ysum[nz - 2])
        / (2.0 * ysum[nz - 1])
        - 1.0
        / (dzi[-1])
        * Dzz[-1]
        / 2.0
        * (
            -1.0 / Hpi[-1]
            + ms * g[-1] / (Navo * kb * Ti[-1])
            + alpha / Ti[-1] * (Tco[-1] - Tco[-2]) / dzi[-1]
        )
        + ((vs[-1] < 0) * vs[-1]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[(nz - 1) - 1]] -= (
        1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[(nz - 1) - 1] + ysum[nz - 1])
        / (2.0 * ysum[(nz - 1) - 1])
        + ((vz[-1] > 0) * vz[-1]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[(nz - 1) - 1]] -= (
        1.0
        / (dzi[nz - 2])
        * (Dzz[nz - 2] / dzi[nz - 2])
        * (ysum[nz - 1] + ysum[nz - 2])
        / (2.0 * ysum[(nz - 1) - 1])
        - 1.0
        / (dzi[-1])
        * Dzz[-1]
        / 2.0
        * (
            -1.0 / Hpi[-1]
            + ms * g[-1] / (Navo * kb * Ti[-1])
            + alpha / Ti[-1] * (Tco[-1] - Tco[-2]) / dzi[-1]
        )
        + ((vs[-1] > 0) * vs[-1]) / dzi[-1]
    )

    return dfdy

`lhs_jac_settling_vm(var, atm)`

directly constructing lhs = 1./(rh)sparse.identity(ni*nz) - dfdy jacobian matrix for dn/dt + dphi/dz = P - L (including molecular diffusion and gravitation settling for particles) zero-flux BC: 1st derivitive of y is zero inc. vs from molecular diffusion

Source code in src/vulcan/op.py

def lhs_jac_settling_vm(self, var, atm):
    """
    directly constructing lhs = 1./(r*h)*sparse.identity(ni*nz) - dfdy
    jacobian matrix for dn/dt + dphi/dz = P - L (including molecular diffusion and gravitation settling for particles)
    zero-flux BC:  1st derivitive of y is zero
    inc. vs from molecular diffusion
    """
    y = var.y.copy()
    nz = self.cfg.nz

    # TEST condensation excluding non-gaseous species
    if self.cfg.non_gas_sp:
        ysum = np.sum(y[:, atm.gas_indx], axis=1)
    else:
        ysum = np.sum(y, axis=1)
    # TEST condensation excluding non-gaseous species
    dzi = atm.dzi.copy()
    Kzz = atm.Kzz.copy()
    Dzz = atm.Dzz.copy()
    vz = atm.vz.copy()
    vs = atm.vs.copy()
    alpha = atm.alpha.copy()
    Tco = atm.Tco.copy()
    mu, ms = atm.mu.copy(), atm.ms.copy()
    g = atm.g
    vm = atm.vm

    Ti = atm.Ti.copy()
    Hpi = atm.Hpi.copy()

    # c0 = 1./(r*h) where r = 1. + 1./2.**0.5
    r = 1.0 + 1.0 / 2.0**0.5
    c0 = 1.0 / (r * var.dt)
    dfdy = neg_achemjac(y, atm.M, var.k)
    np.fill_diagonal(dfdy, c0 + np.diag(dfdy))
    j_indx = []

    for j in range(nz):
        j_indx.append(np.arange(j * ni, j * ni + ni))

    for j in range(1, nz - 1):
        # excluding the buttom and the top cell
        # at j level consists of ni species
        dz_ave = 0.5 * (dzi[j - 1] + dzi[j])
        dfdy[j_indx[j], j_indx[j]] -= (
            -1.0
            / dz_ave
            * (
                Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Kzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / 2.0
            )
            / ysum[j]
            - ((vz[j] > 0) * vz[j] - (vz[j - 1] < 0) * vz[j - 1]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j + 1]] -= (
            1.0 / dz_ave * (Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / (2.0 * ysum[j + 1]))
            - ((vz[j] < 0) * vz[j]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j - 1]] -= (
            1.0
            / dz_ave
            * (Kzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / (2.0 * ysum[j - 1]))
            + ((vz[j - 1] > 0) * vz[j - 1]) / dz_ave
        )

        # [j_indx[j], j_indx[j]] has size ni*ni
        dfdy[j_indx[j], j_indx[j]] -= (
            -1.0
            / dz_ave
            * (
                Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Dzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / 2.0
            )
            / ysum[j]
            - ((vs[j] > 0) * vs[j] - (vs[j - 1] < 0) * vs[j - 1]) / dz_ave
            - ((vm[j] > 0) * vm[j] - (vm[j - 1] < 0) * vm[j - 1]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j + 1]] -= (
            1.0 / dz_ave * (Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / (2.0 * ysum[j + 1]))
            - ((vs[j] < 0) * vs[j]) / dz_ave
            - ((vm[j] < 0) * vm[j]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j - 1]] -= (
            1.0
            / dz_ave
            * (Dzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / (2.0 * ysum[j - 1]))
            + ((vs[j - 1] > 0) * vs[j - 1]) / dz_ave
            + ((vm[j - 1] > 0) * vm[j - 1]) / dz_ave
        )

    dfdy[j_indx[0], j_indx[0]] -= (
        -1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[0])
        - ((vz[0] > 0) * vz[0]) / dzi[0]
    )
    dfdy[j_indx[0], j_indx[0]] -= (
        -1.0 / (dzi[0]) * (Dzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[0])
        - ((vs[0] > 0) * vs[0]) / dzi[0]
    )
    # deposition velocity
    if self.cfg.use_botflux:
        dfdy[j_indx[0], j_indx[0]] -= -1.0 * atm.bot_vdep / dzi[0]

    # diffusion-limited escape
    if self.cfg.diff_esc:  # not empty list
        diff_lim = np.zeros(ni)
        for sp in self.cfg.diff_esc:
            if y[-1, species.index(sp)] > 0:
                diff_lim[species.index(sp)] += (
                    atm.top_flux[species.index(sp)] / y[-1, species.index(sp)]
                )
        dfdy[j_indx[-1], j_indx[-1]] -= diff_lim  # negative

    dfdy[j_indx[0], j_indx[1]] -= (
        1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[1])
        - ((vz[0] < 0) * vz[0]) / dzi[0]
    )
    dfdy[j_indx[0], j_indx[1]] -= (
        1.0 / (dzi[0]) * (Dzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[1])
        - ((vs[0] < 0) * vs[0]) / dzi[0]
    )

    dfdy[j_indx[nz - 1], j_indx[nz - 1]] -= (
        -1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[(nz - 1) - 1] + ysum[nz - 1])
        / (2.0 * ysum[nz - 1])
        + ((vz[-1] < 0) * vz[-1]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[nz - 1]] -= (
        -1.0
        / (dzi[nz - 2])
        * (Dzz[nz - 2] / dzi[nz - 2])
        * (ysum[nz - 1] + ysum[nz - 2])
        / (2.0 * ysum[nz - 1])
        + ((vs[-1] < 0) * vs[-1]) / dzi[-1]
        + ((vm[-1] < 0) * vm[-1]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[(nz - 1) - 1]] -= (
        1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[(nz - 1) - 1] + ysum[nz - 1])
        / (2.0 * ysum[(nz - 1) - 1])
        + ((vz[-1] > 0) * vz[-1]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[(nz - 1) - 1]] -= (
        1.0
        / (dzi[nz - 2])
        * (Dzz[nz - 2] / dzi[nz - 2])
        * (ysum[nz - 1] + ysum[nz - 2])
        / (2.0 * ysum[(nz - 1) - 1])
        + ((vs[-1] > 0) * vs[-1]) / dzi[-1]
        + ((vm[-1] > 0) * vm[-1]) / dzi[-1]
    )

    return dfdy

`lhs_jac_tot(var, atm)`

directly constructing lhs = 1./(rh)sparse.identity(ni*nz) - dfdy jacobian matrix for dn/dt + dphi/dz = P - L (including molecular diffusion) zero-flux BC: 1st derivitive of y is zero

Source code in src/vulcan/op.py

def lhs_jac_tot(self, var, atm):
    """
    directly constructing lhs = 1./(r*h)*sparse.identity(ni*nz) - dfdy
    jacobian matrix for dn/dt + dphi/dz = P - L (including molecular diffusion)
    zero-flux BC:  1st derivitive of y is zero
    """

    nz = self.cfg.nz

    y = var.y.copy()
    # TEST condensation excluding non-gaseous species
    if self.cfg.use_condense:
        ysum = np.sum(y[:, atm.gas_indx], axis=1)
        # ysum = np.sum(y, axis=1)
    else:
        ysum = np.sum(y, axis=1)
    # TEST condensation excluding non-gaseous species
    dzi = atm.dzi.copy()
    Kzz = atm.Kzz.copy()
    Dzz = atm.Dzz.copy()
    vz = atm.vz.copy()
    alpha = atm.alpha.copy()
    Tco = atm.Tco
    mu, ms = atm.mu.copy(), atm.ms.copy()
    g = atm.g

    Ti = atm.Ti.copy()
    Hpi = atm.Hpi.copy()

    # c0 = 1./(r*h) where r = 1. + 1./2.**0.5
    r = 1.0 + 1.0 / 2.0**0.5
    c0 = 1.0 / (r * var.dt)
    dfdy = neg_achemjac(y, atm.M, var.k, self.cfg.nz)
    np.fill_diagonal(dfdy, c0 + np.diag(dfdy))
    j_indx = []

    for j in range(nz):
        j_indx.append(np.arange(j * ni, j * ni + ni))

    for j in range(1, nz - 1):
        # excluding the buttom and the top cell
        # at j level consists of ni species
        dz_ave = 0.5 * (dzi[j - 1] + dzi[j])
        dfdy[j_indx[j], j_indx[j]] -= (
            -1.0
            / dz_ave
            * (
                Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Kzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / 2.0
            )
            / ysum[j]
            - ((vz[j] > 0) * vz[j] - (vz[j - 1] < 0) * vz[j - 1]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j + 1]] -= (
            1.0 / dz_ave * (Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / (2.0 * ysum[j + 1]))
            - ((vz[j] < 0) * vz[j]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j - 1]] -= (
            1.0
            / dz_ave
            * (Kzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / (2.0 * ysum[j - 1]))
            + ((vz[j - 1] > 0) * vz[j - 1]) / dz_ave
        )

        # [j_indx[j], j_indx[j]] has size ni*ni
        dfdy[j_indx[j], j_indx[j]] -= -1.0 / dz_ave * (
            Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
            + Dzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / 2.0
        ) / ysum[j] + 1.0 / (2.0 * dz_ave) * (
            Dzz[j]
            * (
                -1.0 / Hpi[j]
                + ms * g[j] / (Navo * kb * Ti[j])
                + alpha / Ti[j] * (Tco[j + 1] - Tco[j]) / dzi[j]
            )
            - Dzz[j - 1]
            * (
                -1.0 / Hpi[j - 1]
                + ms * g[j] / (Navo * kb * Ti[j - 1])
                + alpha / Ti[j - 1] * (Tco[j] - Tco[j - 1]) / dzi[j - 1]
            )
        )
        dfdy[j_indx[j], j_indx[j + 1]] -= 1.0 / dz_ave * (
            Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / (2.0 * ysum[j + 1])
        ) + 1.0 / (2.0 * dz_ave) * Dzz[j] * (
            -1.0 / Hpi[j]
            + ms * g[j + 1] / (Navo * kb * Ti[j])
            + alpha / Ti[j] * (Tco[j + 1] - Tco[j]) / dzi[j]
        )
        dfdy[j_indx[j], j_indx[j - 1]] -= 1.0 / dz_ave * (
            Dzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / (2.0 * ysum[j - 1])
        ) - 1.0 / (2.0 * dz_ave) * Dzz[j - 1] * (
            -1.0 / Hpi[j - 1]
            + ms * g[j - 1] / (Navo * kb * Ti[j - 1])
            + alpha / Ti[j - 1] * (Tco[j] - Tco[j - 1]) / dzi[j - 1]
        )

    dfdy[j_indx[0], j_indx[0]] -= (
        -1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[0])
        - ((vz[0] > 0) * vz[0]) / dzi[0]
    )
    dfdy[j_indx[0], j_indx[0]] -= -1.0 / (dzi[0]) * (Dzz[0] / dzi[0]) * (
        ysum[1] + ysum[0]
    ) / (2.0 * ysum[0]) + 1.0 / (dzi[0]) * Dzz[0] / 2.0 * (
        -1.0 / Hpi[0]
        + ms * g[0] / (Navo * kb * Ti[0])
        + alpha / Ti[0] * (Tco[1] - Tco[0]) / dzi[0]
    )
    # deposition velocity
    if self.cfg.use_botflux:
        dfdy[j_indx[0], j_indx[0]] -= -1.0 * atm.bot_vdep / dzi[0]

    dfdy[j_indx[0], j_indx[1]] -= (
        1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[1])
        - ((vz[0] < 0) * vz[0]) / dzi[0]
    )
    dfdy[j_indx[0], j_indx[1]] -= 1.0 / (dzi[0]) * (Dzz[0] / dzi[0]) * (
        ysum[1] + ysum[0]
    ) / (2.0 * ysum[1]) + 1.0 / (dzi[0]) * Dzz[0] / 2.0 * (
        -1.0 / Hpi[0]
        + ms * g[0] / (Navo * kb * Ti[0])
        + alpha / Ti[0] * (Tco[1] - Tco[0]) / dzi[0]
    )

    dfdy[j_indx[nz - 1], j_indx[nz - 1]] -= (
        -1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[(nz - 1) - 1] + ysum[nz - 1])
        / (2.0 * ysum[nz - 1])
        + ((vz[-1] < 0) * vz[-1]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[nz - 1]] -= -1.0 / (dzi[nz - 2]) * (
        Dzz[nz - 2] / dzi[nz - 2]
    ) * (ysum[nz - 1] + ysum[nz - 2]) / (2.0 * ysum[nz - 1]) - 1.0 / (dzi[-1]) * Dzz[
        -1
    ] / 2.0 * (
        -1.0 / Hpi[-1]
        + ms * g[-1] / (Navo * kb * Ti[-1])
        + alpha / Ti[-1] * (Tco[-1] - Tco[-2]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[(nz - 1) - 1]] -= (
        1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[(nz - 1) - 1] + ysum[nz - 1])
        / (2.0 * ysum[(nz - 1) - 1])
        + ((vz[-1] > 0) * vz[-1]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[(nz - 1) - 1]] -= 1.0 / (dzi[nz - 2]) * (
        Dzz[nz - 2] / dzi[nz - 2]
    ) * (ysum[nz - 1] + ysum[nz - 2]) / (2.0 * ysum[(nz - 1) - 1]) - 1.0 / (dzi[-1]) * Dzz[
        -1
    ] / 2.0 * (
        -1.0 / Hpi[-1]
        + ms * g[-1] / (Navo * kb * Ti[-1])
        + alpha / Ti[-1] * (Tco[-1] - Tco[-2]) / dzi[-1]
    )

    return dfdy

`lhs_jac_tot_vm(var, atm)`

directly constructing lhs = 1./(rh)sparse.identity(ni*nz) - dfdy jacobian matrix for dn/dt + dphi/dz = P - L (including molecular diffusion) zero-flux BC: 1st derivitive of y is zero inc. vm from molecular diffusion

Source code in src/vulcan/op.py

def lhs_jac_tot_vm(self, var, atm):
    """
    directly constructing lhs = 1./(r*h)*sparse.identity(ni*nz) - dfdy
    jacobian matrix for dn/dt + dphi/dz = P - L (including molecular diffusion)
    zero-flux BC:  1st derivitive of y is zero
    inc. vm from molecular diffusion
    """
    y = var.y.copy()
    nz = self.cfg.nz

    # TEST condensation excluding non-gaseous species
    if self.cfg.use_condense:
        ysum = np.sum(y[:, atm.gas_indx], axis=1)
        # ysum = np.sum(y, axis=1)
    else:
        ysum = np.sum(y, axis=1)
    # TEST condensation excluding non-gaseous species
    dzi = atm.dzi.copy()
    Kzz = atm.Kzz.copy()
    Dzz = atm.Dzz.copy()
    vz = atm.vz.copy()
    alpha = atm.alpha.copy()
    Tco = atm.Tco.copy()
    mu, ms = atm.mu.copy(), atm.ms.copy()
    g = atm.g
    vm = atm.vm

    Ti = atm.Ti.copy()
    Hpi = atm.Hpi.copy()

    # c0 = 1./(r*h) where r = 1. + 1./2.**0.5
    r = 1.0 + 1.0 / 2.0**0.5
    c0 = 1.0 / (r * var.dt)
    dfdy = neg_achemjac(y, atm.M, var.k)
    np.fill_diagonal(dfdy, c0 + np.diag(dfdy))
    j_indx = []

    for j in range(nz):
        j_indx.append(np.arange(j * ni, j * ni + ni))

    for j in range(1, nz - 1):
        # excluding the buttom and the top cell
        # at j level consists of ni species
        dz_ave = 0.5 * (dzi[j - 1] + dzi[j])
        dfdy[j_indx[j], j_indx[j]] -= (
            -1.0
            / dz_ave
            * (
                Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Kzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / 2.0
            )
            / ysum[j]
            - ((vz[j] > 0) * vz[j] - (vz[j - 1] < 0) * vz[j - 1]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j + 1]] -= (
            1.0 / dz_ave * (Kzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / (2.0 * ysum[j + 1]))
            - ((vz[j] < 0) * vz[j]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j - 1]] -= (
            1.0
            / dz_ave
            * (Kzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / (2.0 * ysum[j - 1]))
            + ((vz[j - 1] > 0) * vz[j - 1]) / dz_ave
        )

        # [j_indx[j], j_indx[j]] has size ni*ni
        dfdy[j_indx[j], j_indx[j]] -= (
            -1.0
            / dz_ave
            * (
                Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / 2.0
                + Dzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / 2.0
            )
            / ysum[j]
            - ((vm[j] > 0) * vm[j] - (vm[j - 1] < 0) * vm[j - 1]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j + 1]] -= (
            1.0 / dz_ave * (Dzz[j] / dzi[j] * (ysum[j + 1] + ysum[j]) / (2.0 * ysum[j + 1]))
            - ((vm[j] < 0) * vm[j]) / dz_ave
        )
        dfdy[j_indx[j], j_indx[j - 1]] -= (
            1.0
            / dz_ave
            * (Dzz[j - 1] / dzi[j - 1] * (ysum[j - 1] + ysum[j]) / (2.0 * ysum[j - 1]))
            + ((vm[j - 1] > 0) * vm[j - 1]) / dz_ave
        )

    dfdy[j_indx[0], j_indx[0]] -= (
        -1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[0])
        - ((vz[0] > 0) * vz[0]) / dzi[0]
    )
    dfdy[j_indx[0], j_indx[0]] -= (
        -1.0 / (dzi[0]) * (Dzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[0])
        - ((vm[0] > 0) * vm[0]) / dzi[0]
    )
    # deposition velocity
    if self.cfg.use_botflux:
        dfdy[j_indx[0], j_indx[0]] -= -1.0 * atm.bot_vdep / dzi[0]
    # diffusion-limited escape
    if self.cfg.diff_esc:  # not empty list
        diff_lim = np.zeros(ni)
        for sp in self.cfg.diff_esc:
            if y[-1, species.index(sp)] > 0:
                diff_lim[species.index(sp)] += (
                    atm.top_flux[species.index(sp)] / y[-1, species.index(sp)]
                )
        dfdy[j_indx[-1], j_indx[-1]] -= diff_lim  # negative

    dfdy[j_indx[0], j_indx[1]] -= (
        1.0 / (dzi[0]) * (Kzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[1])
        - ((vz[0] < 0) * vz[0]) / dzi[0]
    )
    dfdy[j_indx[0], j_indx[1]] -= (
        1.0 / (dzi[0]) * (Dzz[0] / dzi[0]) * (ysum[1] + ysum[0]) / (2.0 * ysum[1])
        - ((vm[0] < 0) * vm[0]) / dzi[0]
    )

    dfdy[j_indx[nz - 1], j_indx[nz - 1]] -= (
        -1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[(nz - 1) - 1] + ysum[nz - 1])
        / (2.0 * ysum[nz - 1])
        + ((vz[-1] < 0) * vz[-1]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[nz - 1]] -= (
        -1.0
        / (dzi[nz - 2])
        * (Dzz[nz - 2] / dzi[nz - 2])
        * (ysum[nz - 1] + ysum[nz - 2])
        / (2.0 * ysum[nz - 1])
        + ((vm[-1] < 0) * vm[-1]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[(nz - 1) - 1]] -= (
        1.0
        / (dzi[nz - 2])
        * (Kzz[nz - 2] / dzi[nz - 2])
        * (ysum[(nz - 1) - 1] + ysum[nz - 1])
        / (2.0 * ysum[(nz - 1) - 1])
        + ((vz[-1] > 0) * vz[-1]) / dzi[-1]
    )
    dfdy[j_indx[nz - 1], j_indx[(nz - 1) - 1]] -= (
        1.0
        / (dzi[nz - 2])
        * (Dzz[nz - 2] / dzi[nz - 2])
        * (ysum[nz - 1] + ysum[nz - 2])
        / (2.0 * ysum[(nz - 1) - 1])
        + ((vm[-1] > 0) * vm[-1]) / dzi[-1]
    )

    return dfdy

`loss(data_var)`

Source code in src/vulcan/op.py

def loss(self, data_var):

    y = data_var.y
    atom_list = self.cfg.atom_list

    # changed atom_tot to dictionary atom_sum
    atom_sum = data_var.atom_sum

    for atom in atom_list:
        # data_var.atom_sum[atom] = np.sum([compo[compo_row.index(species[i])][atom] * data_var.y[:,i] for i in range(ni)])
        # TEST V scaling
        if atom not in getattr(self.cfg, 'loss_ex', []):  # shami added 2024
            data_var.atom_sum[atom] = np.sum(
                [
                    compo[compo_row.index(species[i])][atom] * data_var.y[:, i]
                    for i in range(ni)
                ]
            )  # *data_var.v_ratio
            data_var.atom_loss[atom] = (
                data_var.atom_sum[atom] - data_var.atom_ini[atom]
            ) / data_var.atom_ini[atom]

    return data_var

`print_lossBig(para)`

Source code in src/vulcan/op.py

def print_lossBig(self, para):

    log.warning('Element conservation is violated too large')
    log.warning('at step: ' + str(para.count))
    log.warning(spacer)

`print_nega(data_var, data_para)`

Source code in src/vulcan/op.py

def print_nega(self, data_var, data_para):

    nega_i = np.where(data_var.y < 0)
    log.warning(
        'Negative y at time '
        + str('{:.2e}'.format(data_var.t))
        + ' and step: '
        + str(data_para.count)
    )
    log.warning('Negative values:' + str(data_var.y[data_var.y < 0]))
    log.warning('from levels: ' + str(nega_i[0]))
    log.warning('species: ' + str([species[s] for s in nega_i[1]]))
    log.warning('dt= ' + str(data_var.dt))
    log.warning('...reset dt to dt*0.2...')
    log.warning(spacer)

`reset_y(var, dt_reduc=0.5)`

reset y and reduce dt by dt_reduc

Source code in src/vulcan/op.py

def reset_y(self, var, dt_reduc=0.5):
    """
    reset y and reduce dt by dt_reduc
    """

    # reset and store y and dt
    var.y = var.y_prev
    var.dt *= dt_reduc
    # var.dt = np.maximum(var.dt, self.cfg.dt_min)

    return var

`step_ok(var, para, loss_eps=0.1, rtol=0.6)`

Source code in src/vulcan/op.py

def step_ok(self, var, para, loss_eps=1e-1, rtol=0.6):
    if (
        np.all(var.y >= 0)
        and np.amax(
            np.abs(
                np.fromiter(var.atom_loss.values(), float)
                - np.fromiter(var.atom_loss_prev.values(), float)
            )
        )
        < loss_eps
        and para.delta <= rtol
    ):
        return True
    else:
        return False

`step_reject(var, para, loss_eps=0.1, rtol=0.6)`

Source code in src/vulcan/op.py

def step_reject(self, var, para, loss_eps=1e-1, rtol=0.6):

    if para.delta > rtol:  # truncation error larger than the tolerence value
        para.delta_count += 1

    elif np.any(var.y < 0):
        para.nega_count += 1
        if self.cfg.use_print_prog:
            self.print_nega(
                var, para
            )  # print the info for the negative solutions (where y < 0)
        # print input: y, t, count, dt

    else:  # meaning np.amax( np.abs( np.abs(y_loss) - np.abs(loss_prev) ) )<loss_eps
        para.loss_count += 1
        if self.cfg.use_print_prog:
            self.print_lossBig(para)

    # reset y and dt to the values at previous step
    var = self.reset_y(var, dt_reduc=self.cfg.dt_var_min)

    if var.dt < self.cfg.dt_min:
        var.dt = self.cfg.dt_min
        var.y[var.y < 0] = 0.0  # clipping of negative values

        log.warning(
            'Keep producing negative values! Clipping negative solutions and moving on!'
        )
        return True

    return False

`thomas_vec(a, b, c, d)`

Thomas vectorized solver, a b c d refer to http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm d is a matrix not used in this current version

Source code in src/vulcan/op.py

def thomas_vec(a, b, c, d):
    """
    Thomas vectorized solver, a b c d refer to http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
    d is a matrix
    not used in this current version
    """
    # number of equations
    nf = len(a)
    aa, bb, cc, dd = map(np.copy, (a, b, c, d))
    # d needs to reshape
    dd = dd.reshape(nf, -1)
    # C' and D'
    cp = [cc[0] / bb[0]]
    dp = [dd[0] / bb[0]]
    x = np.zeros((nf, np.shape(dd)[1]))

    for i in range(1, nf - 1):
        cp.append(cc[i] / (bb[i] - aa[i] * cp[i - 1]))
        dp.append((dd[i] - aa[i] * dp[i - 1]) / (bb[i] - aa[i] * cp[i - 1]))

    dp.append(
        (dd[(nf - 1)] - aa[(nf - 1)] * dp[(nf - 1) - 1])
        / (bb[(nf - 1)] - aa[(nf - 1)] * cp[(nf - 1) - 1])
    )  # nf-1 is the last element
    x[nf - 1] = dp[nf - 1] / 1
    for i in range(nf - 2, -1, -1):
        x[i] = dp[i] - cp[i] * x[i + 1]

    return x

`Ros2(vulcan_cfg)`

Bases: ODESolver

class inheritance from ODEsolver for 2nd order Rosenbrock solver

Source code in src/vulcan/op.py

def __init__(self, vulcan_cfg: Config):
    # ODESolver.__init__(self)
    super().__init__(vulcan_cfg)

`naming_solver(para)`

Source code in src/vulcan/op.py

def naming_solver(self, para):
    if self.cfg.use_moldiff:
        log.info('Include molecular diffusion.')
    else:
        log.info('No molecular diffusion.')
    para.solver_str = 'solver'

`one_step(var, atm, para)`

Source code in src/vulcan/op.py

def one_step(self, var, atm, para):

    while True:
        var, para = getattr(self, para.solver_str)(var, atm, para)

        # clipping small negative values and also calculating atomic loss (atom_loss)
        var, para = self.clip(
            var,
            para,
            atm,
            pos_cut=self.cfg.pos_cut,
            nega_cut=self.cfg.nega_cut,
        )

        if self.step_ok(var, para, loss_eps=self.cfg.loss_eps, rtol=self.cfg.rtol):
            break
        elif self.step_reject(var, para, loss_eps=self.cfg.loss_eps, rtol=self.cfg.rtol):
            break

    return var, para

`solver(var, atm, para)`

2nd order Rosenbrock [Verwer et al. 1997] with banded-matrix solver with switches to include the molecular diffusion or not

Source code in src/vulcan/op.py

def solver(self, var, atm, para):
    """
    2nd order Rosenbrock [Verwer et al. 1997] with banded-matrix solver
    with switches to include the molecular diffusion or not
    """

    nz = self.cfg.nz
    y, ymix, h, k = var.y, var.ymix, var.dt, var.k
    M, dzi, Kzz = atm.M, atm.dzi, atm.Kzz

    if not self.cfg.use_vm_mol:
        if self.cfg.use_moldiff and not self.cfg.use_settling:
            diffdf = self.diffdf
            jac_tot = self.lhs_jac_tot
        elif self.cfg.use_moldiff and self.cfg.use_settling:
            diffdf = self.diffdf_settling
            jac_tot = self.lhs_jac_settling
        else:
            diffdf = self.diffdf_no_mol
            jac_tot = self.lhs_jac_no_mol
    else:  # vulcan_cfg.use_vm_mol :
        if self.cfg.use_moldiff and not self.cfg.use_settling:
            diffdf = self.diffdf_vm
            jac_tot = self.lhs_jac_tot_vm
        elif self.cfg.use_moldiff and self.cfg.use_settling:
            diffdf = self.diffdf_settling_vm
            jac_tot = self.lhs_jac_settling_vm
        else:
            diffdf = self.diffdf_no_mol
            jac_tot = self.lhs_jac_no_mol

    r = 1.0 + 1.0 / 2.0**0.5

    df = chemdf(y, M, k).flatten() + diffdf(y, atm).flatten()
    lhs = jac_tot(var, atm)

    # Fixed species including only below the cold trap # TEST 2022
    if self.cfg.use_condense and para.fix_species_start:
        for sp in self.cfg.fix_species:
            if (
                not self.cfg.fix_species_from_coldtrap_lev
            ):  # if Ptop is not specified, fix the whole column # TEST2022
                pass
            else:
                pfix_indx = atm.conden_min_lev[sp]
                atm.fix_sp_indx[sp] = np.arange(
                    species.index(sp), species.index(sp) + ni * (pfix_indx), ni
                )

            df[atm.fix_sp_indx[sp]] = 0
            lhs[atm.fix_sp_indx[sp], :] = 0
            lhs[atm.fix_sp_indx[sp], atm.fix_sp_indx[sp]] = (
                1.0 / (r * h)
            )  # cuz the jacobian func is directly outputing 1./(r*h)*sparse.identity(ni*nz) - dfdy

    if self.cfg.use_ion:
        df[atm.fix_e_indx] = 0
        lhs[atm.fix_e_indx, :] = 0
        lhs[atm.fix_e_indx, atm.fix_e_indx] = 1.0 / (r * h)

    lhs_b, bw = self.store_bandM(lhs, ni, nz)
    k1_flat = scipy.linalg.solve_banded((bw, bw), lhs_b, df)
    k1 = k1_flat.reshape(y.shape)

    yk2 = y + k1 / r
    df = chemdf(yk2, M, k).flatten() + diffdf(yk2, atm).flatten()

    # TEST condensation
    # Fixed species
    if self.cfg.use_condense and para.fix_species_start:
        for sp in self.cfg.fix_species:
            df[atm.fix_sp_indx[sp]] = 0
    if self.cfg.use_ion:
        df[atm.fix_e_indx] = 0

    rhs = df - 2.0 / (r * h) * k1_flat
    k2 = scipy.linalg.solve_banded((bw, bw), lhs_b, rhs)
    k2 = k2.reshape(y.shape)

    sol = y + 3.0 / (2.0 * r) * k1 + 1 / (2.0 * r) * k2

    ### for Hycean ###
    if (
        getattr(self.cfg, 'use_fix_H2He', False)
        and 'H2' not in self.cfg.use_fix_sp_bot
        and var.t > 1e6
    ):
        self.cfg.use_fix_sp_bot['H2'] = var.ymix[0, species.index('H2')]
        self.cfg.use_fix_sp_bot['He'] = var.ymix[0, species.index('He')]
        print(
            'After 1e6 sec, H2 and He are fixed at '
            + str((var.ymix[0, species.index('H2')], var.ymix[0, species.index('He')]))
        )

        self.fix_sp_bot_index = [species.index(sp) for sp in self.cfg.use_fix_sp_bot.keys()]
        self.fix_sp_bot_mix = np.array(
            [self.cfg.use_fix_sp_bot[sp] for sp in self.cfg.use_fix_sp_bot.keys()]
        )
    ### for Hycean ###

    # setting particles on the surace = 0
    if self.cfg.use_fix_sp_bot:  # if use_fix_sp_bot = {} (empty), it returns false
        sol[0, self.fix_sp_bot_index] = self.fix_sp_bot_mix * atm.n_0[0]

    delta = np.abs(sol - yk2)
    delta[ymix < self.mtol] = 0
    delta[sol < self.atol] = 0

    # neglecting the errors at the surface
    if self.cfg.use_botflux or self.cfg.use_fix_sp_bot:
        delta[0] = 0

    # TEST condensation 2022
    if self.cfg.use_condense:
        delta[:, self.non_gas_sp_index] = 0
        delta[:, self.condense_sp_index] = 0

        if para.fix_species_start:
            for sp in self.cfg.fix_species:
                if (
                    not self.cfg.fix_species_from_coldtrap_lev
                ):  # if Ptop is not specified, fix the whole column # TEST2022
                    sol[:, species.index(sp)] = var.fix_y[sp].copy()
                else:
                    # pfix_indx = min( range(len(atm.pco)), key=lambda i: abs(atm.pco[i]- self.cfg.fix_species_Ptop[0] ))
                    pfix_indx = atm.conden_min_lev[sp]
                    sol[:pfix_indx, species.index(sp)] = var.fix_y[sp].copy()[:pfix_indx]

                delta[:, species.index(sp)] = 0

    if self.cfg.use_print_delta and para.count % self.cfg.print_prog_num == 0:
        max_indx = np.nanargmax(delta / sol, axis=1)
        max_lev_indx = np.nanargmax(delta / sol)
        log.info(
            'Largest delta (truncation error) from nz = ' + str(int(max_lev_indx / ni))
        )
        log.info(np.array(species)[max_indx])
        log.info(
            'Largest delta (truncation error) from '
            + species[max_indx % ni]
            + ' at nz = '
            + str(int(max_indx / ni))
        )

    delta = np.amax(delta[sol > 0] / sol[sol > 0])

    var.y = sol

    # # TEST condensation excluding non-gaseous species
    if self.cfg.non_gas_sp:
        var.ymix = var.y / np.vstack(np.sum(var.y[:, atm.gas_indx], axis=1))
    else:
        var.ymix = var.y / np.vstack(np.sum(var.y, axis=1))
    # TEST condensation excluding non-gaseous species

    para.delta = delta

    # use charge balance to obtain the number density of electrons (such that [ions] = [e])
    if self.cfg.use_ion:
        # clear e
        var.y[:, species.index('e')] = 0
        # set e such that the net chare is zero
        for sp in var.charge_list:
            var.y[:, species.index('e')] -= (
                compo[compo_row.index(sp)]['e'] * var.y[:, species.index(sp)]
            )

    return var, para

`solver_fix_all_bot(var, atm, para)`

2nd order Rosenbrock [Verwer et al. 1997] with banded-matrix solver with switches to include the molecular diffusion or not

Source code in src/vulcan/op.py

def solver_fix_all_bot(self, var, atm, para):
    """
    2nd order Rosenbrock [Verwer et al. 1997] with banded-matrix solver
    with switches to include the molecular diffusion or not
    """

    nz = self.cfg.nz

    y, ymix, h, k = var.y, var.ymix, var.dt, var.k
    M, dzi, Kzz = atm.M, atm.dzi, atm.Kzz

    # store the fixed bottom level
    bottom = np.copy(ymix[0])

    if self.cfg.use_moldiff:
        diffdf = self.diffdf
        jac_tot = self.lhs_jac_fix_all_bot
    else:
        diffdf = self.diffdf_no_mol
        jac_tot = self.lhs_jac_no_mol_fix_all_bot

    r = 1.0 + 1.0 / 2.0**0.5

    df = chemdf(y, M, k).flatten() + diffdf(y, atm).flatten()
    lhs = jac_tot(var, atm)

    lhs_b, bw = self.store_bandM(lhs, ni, nz)
    k1_flat = scipy.linalg.solve_banded((bw, bw), lhs_b, df)

    k1 = k1_flat.reshape(y.shape)

    yk2 = y + k1 / r
    df = chemdf(yk2, M, k).flatten() + diffdf(yk2, atm).flatten()

    rhs = df - 2.0 / (r * h) * k1_flat
    k2 = scipy.linalg.solve_banded((bw, bw), lhs_b, rhs)
    k2 = k2.reshape(y.shape)

    sol = y + 3.0 / (2.0 * r) * k1 + 1 / (2.0 * r) * k2

    # fixed the bottom layer to yini (in chemical EQ)
    sol[0] = bottom * atm.n_0[0]

    delta = np.abs(sol - yk2)
    delta[ymix < self.mtol] = 0
    delta[sol < self.atol] = 0

    delta = np.amax(delta[sol > 0] / sol[sol > 0])

    var.y = sol

    # # TEST condensation excluding non-gaseous species
    if self.cfg.non_gas_sp:
        var.ymix = var.y / np.vstack(np.sum(var.y[:, atm.gas_indx], axis=1))
    else:
        var.ymix = var.y / np.vstack(np.sum(var.y, axis=1))
    # TEST condensation excluding non-gaseous species

    para.delta = delta

    # use charge balance to obtain the number density of electrons (such that [ions] = [e])
    if self.cfg.use_ion:
        # clear e
        var.y[:, species.index('e')] = 0
        # set e such that the net chare is zero
        for sp in var.charge_list:
            var.y[:, species.index('e')] -= (
                compo[compo_row.index(sp)]['e'] * var.y[:, species.index(sp)]
            )

    return var, para

`step_size(var, para, dt_var_min=2, dt_var_max=0.5, dt_min=1e-10, dt_max=1e+18)`

step-size control by delta(truncation error) for the Rosenbrock method

Source code in src/vulcan/op.py

def step_size(self, var, para, dt_var_min=2, dt_var_max=0.5, dt_min=1e-10, dt_max=1e18):
    """
    step-size control by delta(truncation error) for the Rosenbrock method
    """
    h = var.dt
    delta = para.delta
    rtol = self.cfg.rtol

    if delta == 0:
        delta = 0.01 * rtol
    h_factor = 0.9 * (rtol / delta) ** 0.5  # 0.9 is simply a safety factor
    h_factor = np.maximum(h_factor, dt_var_min)
    h_factor = np.minimum(h_factor, dt_var_max)

    h *= h_factor
    h = np.maximum(h, dt_min)
    h = np.minimum(h, dt_max)

    # store the adopted dt
    var.dt = h

    return var

`store_bandM(a, nb, nn)`

store block-tridiagonal matrix(bandwidth=1) into diagonal ordered form (http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_banded.html) a : square block-tridiagonal matirx nb: size of the block matrix (number of species) nn: number of the block matrices (number of layers)

Source code in src/vulcan/op.py

def store_bandM(self, a, nb, nn):
    """
    store block-tridiagonal matrix(bandwidth=1) into diagonal ordered form
    (http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_banded.html)
    a : square block-tridiagonal matirx
    nb: size of the block matrix (number of species)
    nn: number of the block matrices (number of layers)
    """

    # band width (treat block-banded as banded matrix)
    bw = 2 * nb - 1
    ab = np.zeros((2 * bw + 1, nb * nn))

    # first 2 columns
    for i in range(0, 2 * nb):
        ab[-(2 * nb + i) :, i] = a[0 : 2 * nb + i, i]

    # middle
    for i in range(2 * nb, nn * nb - 2 * nb):
        ab[:, i] = a[(i - 2 * nb + 1) : (i - 2 * nb + 1) + (2 * bw + 1), i]

    # last 2 columns
    for ne, i in enumerate(range(nn * nb - 2 * nb, nn * nb)):
        ab[: (2 * bw + 1 - ne), i] = a[-(2 * bw + 1 - ne) :, i]

    return (ab, bw)

`Output(vulcan_cfg)`

Bases: object

Source code in src/vulcan/op.py

def __init__(self, vulcan_cfg: Config):

    self.cfg = vulcan_cfg
    output_dir = self.cfg.output_dir + '/'
    out_name = self.cfg.out_name

    if self.cfg.clean_output:
        safe_rm(self.cfg.output_dir)

    if not os.path.exists(output_dir):
        os.makedirs(output_dir)

    if not os.path.exists(self.cfg.plot_dir):
        os.makedirs(self.cfg.plot_dir)

    outfile = output_dir + out_name
    if os.path.isfile(outfile):
        log.warning('Output file already exists. Removing.')
        safe_rm(outfile)

`cfg = vulcan_cfg` `instance-attribute`

`plot_TP(atm)`

Source code in src/vulcan/op.py

def plot_TP(self, atm):
    plot_dir = self.cfg.plot_dir
    # plt.figure('TPK')
    fig, ax1 = plt.subplots()
    ax2 = ax1.twiny()  # ax1 and ax2 share y-axis

    if not self.cfg.plot_height:
        ax1.semilogy(atm.Tco, atm.pco / 1.0e6, c='black')
        ax2.loglog(atm.Kzz, atm.pico[1:-1] / 1.0e6, c='k', ls='--')
        plt.gca().invert_yaxis()
        plt.ylim((self.cfg.P_b / 1.0e6, self.cfg.P_t / 1.0e6))
        ax1.set_ylabel('Pressure (bar)')

    else:  # plotting with height
        ax1.plot(atm.Tco, atm.zmco / 1.0e5, c='black')
        ax2.semilogx(atm.Kzz, atm.zmco[1:] / 1.0e5, c='k', ls='--')
        ax1.set_ylabel('Height (km)')

    # plt.xlabel("Temperature (K)")
    ax1.set_xlabel('Temperature (K)')
    ax2.set_xlabel(r'K$_{zz}$ (cm$^2$s$^{-1}$)')

    fig.savefig(plot_dir + f'TPK_initial.{PLT_FMT}', dpi=self.cfg.plot_dpi)

`plot_evo(var, atm, plot_j=-1, plot_ymin=1e-16, dn=1)`

Plot the evolution of mixing ratios over time.

Arguments

var: variable object containing the time evolution data
atm: atmosphere object containing the atmospheric properties
plot_j: index of the atmospheric layer to plot (default: -1)
plot_ymin: minimum y-axis value for the plot (default: 1e-16)
dn: data point interval for plotting (default: 1, meaning plot all points)

Source code in src/vulcan/op.py

def plot_evo(self, var, atm, plot_j=-1, plot_ymin=1e-16, dn=1):
    """
    Plot the evolution of mixing ratios over time.

    Arguments
    ------------
    - var: variable object containing the time evolution data
    - atm: atmosphere object containing the atmospheric properties
    - plot_j: index of the atmospheric layer to plot (default: -1)
    - plot_ymin: minimum y-axis value for the plot (default: 1e-16)
    - dn: data point interval for plotting (default: 1, meaning plot all points)
    """

    plot_spec = self.cfg.plot_spec
    plot_dir = self.cfg.plot_dir
    fig, ax = plt.subplots(1, 1, figsize=(8, 6))

    ymix_time = np.array(var.y_time / atm.n_0[:, np.newaxis])

    for i, sp in enumerate(self.cfg.plot_spec):
        col = gas_cols.get(sp, plt.cm.rainbow(float(i) / len(plot_spec)))
        lbl = tex_labels.get(sp, sp)
        ax.plot(
            var.t_time[::dn], ymix_time[::dn, plot_j, species.index(sp)], c=col, label=lbl
        )

    ax.set_xscale('log')
    ax.set_yscale('log')
    ax.set_xlabel('Time [s]')
    ax.set_ylabel('Volume mixing ratios')
    ax.set_ylim((plot_ymin, 1.2))
    ax.legend(fontsize=10, labelspacing=0.2, loc='upper left', bbox_to_anchor=(1.0, 1.0))
    fig.savefig(plot_dir + f'evo.{PLT_FMT}', dpi=self.cfg.plot_dpi)
    plt.close('all')

`plot_evo_inter(var, atm, plot_j=-1, plot_ymin=1e-16, dn=1)`

plot the evolution when the code is interrupted

Source code in src/vulcan/op.py

def plot_evo_inter(self, var, atm, plot_j=-1, plot_ymin=1e-16, dn=1):
    """
    plot the evolution when the code is interrupted
    """
    var.t_time = np.array(var.t_time)
    ymix_time = np.array(var.y_time / atm.n_0[:, np.newaxis])

    plot_spec = self.cfg.plot_spec
    plot_dir = self.cfg.plot_dir
    plt.figure('evolution')

    for i, sp in enumerate(self.cfg.plot_spec):
        plt.plot(
            var.t_time[::dn],
            ymix_time[::dn, plot_j, species.index(sp)],
            c=plt.cm.rainbow(float(i) / len(plot_spec)),
            label=sp,
        )

    plt.gca().set_xscale('log')
    plt.gca().set_yscale('log')
    plt.xlabel('time')
    plt.ylabel('mixing ratios')
    plt.ylim((plot_ymin, 1.0))
    plt.legend(frameon=0, prop={'size': 14}, loc='best')
    plt.savefig(plot_dir + f'evo.{PLT_FMT}', dpi=self.cfg.plot_dpi)

`plot_flux_update(var, atm, para)`

Source code in src/vulcan/op.py

def plot_flux_update(self, var, atm, para):

    images = []
    plt.ion()

    # fig.add_subplot(121) fig.add_subplot(122)

    (line1,) = plt.plot(np.sum(var.dflux_u, axis=1), atm.pico / 1.0e6, label='up flux')
    (line2,) = plt.plot(
        np.sum(var.dflux_d, axis=1), atm.pico / 1.0e6, label='down flux', ls='--', lw=1.2
    )
    (line3,) = plt.plot(
        np.sum(var.sflux, axis=1), atm.pico / 1.0e6, label='stellar flux', ls=':', lw=1.5
    )

    images.append((line1, line2))

    plt.title(str(para.count) + ' steps and ' + str('{:.2e}'.format(var.t)) + ' s')
    plt.gca().set_xscale('log')
    plt.gca().set_yscale('log')
    plt.gca().invert_yaxis()
    plt.xlim(xmin=1.0e-8)
    plt.ylim((atm.pico[0] / 1.0e6, atm.pico[-1] / 1.0e6))
    plt.legend(frameon=0, prop={'size': 14}, loc=3)
    plt.xlabel('Diffusive flux')
    plt.ylabel('Pressure (bar)')
    plt.show(block=0)
    plt.pause(0.1)
    if self.cfg.use_flux_movie:
        plt.savefig('plot/movie/flux-' + str(para.count) + '.jpg', dpi=self.cfg.plot_dpi)

    plt.close('all')

`plot_update(var, atm, para, fpath=None)`

Source code in src/vulcan/op.py

def plot_update(self, var, atm, para, fpath=None):

    log.debug('Plotting mixing ratios')

    colors = [
        'b',
        'r',
        'c',
        'm',
        'y',
        'k',
        'orange',
        'pink',
        'grey',
        'darkred',
        'darkblue',
        'salmon',
        'chocolate',
        'mediumspringgreen',
        'steelblue',
        'plum',
        'hotpink',
    ]

    fig, ax = plt.subplots(1, 1, figsize=(8, 6))
    color_index = 0
    for sp in self.cfg.plot_spec:
        if sp in tex_labels:
            sp_lab = tex_labels[sp]
        else:
            sp_lab = sp

        if color_index == len(colors):  # when running out of colors
            colors.append(tuple(np.random.rand(3)))

        if sp in gas_cols.keys():
            color = gas_cols[sp]
        else:
            color = colors[color_index]
            color_index += 1

        if not self.cfg.plot_height:
            (line,) = ax.plot(
                var.ymix[:, species.index(sp)], atm.pco / 1.0e6, color=color, label=sp_lab
            )
            if self.cfg.use_condense and sp in self.cfg.condense_sp:
                ax.plot(
                    atm.sat_mix[sp],
                    atm.pco / 1.0e6,
                    color=color,
                    label=sp_lab + ' sat',
                    ls='--',
                )

            ax.set_yscale('log')
            ax.invert_yaxis()
            ax.set_ylabel('Pressure [bar]')
            ax.set_ylim((self.cfg.P_b / 1.0e6, self.cfg.P_t / 1.0e6))

        else:  # plotting with height
            (line,) = ax.plot(
                var.ymix[:, species.index(sp)], atm.zmco / 1.0e5, color=color, label=sp_lab
            )
            if self.cfg.use_condense and sp in self.cfg.condense_sp:
                ax.plot(
                    atm.sat_mix[sp],
                    atm.zco[1:] / 1.0e5,
                    color=color,
                    label=sp_lab + ' sat',
                    ls='--',
                )

            ax.set_ylim((atm.zco[0] / 1e5, atm.zco[-1] / 1e5))
            ax.set_ylabel('Height [km]')

    title = 'i = %5d steps     t = %.2e seconds' % (para.count, var.t)
    title += '\n dy/dt = %.2e' % (var.longdydt)
    ax.set_title(title)
    ax.set_xscale('log')
    ax.set_xlabel('Volume mixing ratio')
    ax.set_xlim(1e-16, 1.2)
    ax.legend(fontsize=10, labelspacing=0.2, loc='upper left', bbox_to_anchor=(1.0, 1.0))

    # temperature profile
    axt = ax.twiny()
    axt.set_xlabel('Temperature [K]')
    if self.cfg.plot_height:
        axt_yarr = atm.zmco / 1.0e5
    else:
        axt_yarr = atm.pco / 1.0e6
    axt.plot(atm.Tco, axt_yarr, color='k', lw=0.5, ls='dotted')

    if fpath is None:
        save_fpath = os.path.join(self.cfg.plot_dir, f'_recent.{PLT_FMT}')
        copy_fpath = os.path.join(self.cfg.plot_dir, f'{para.pic_count:05d}.{PLT_FMT}')
    else:
        save_fpath = fpath
        copy_fpath = None

    log.debug(f'Plotting to {save_fpath}')
    fig.savefig(save_fpath, dpi=self.cfg.plot_dpi, bbox_inches='tight')

    if copy_fpath:
        shutil.copyfile(save_fpath, copy_fpath)

    para.pic_count += 1

`print_end_msg(var, para)`

Source code in src/vulcan/op.py

def print_end_msg(self, var, para):
    log.info('Total atom loss:')
    for atom in self.cfg.atom_list:
        log.info(atom + ': ' + str(var.atom_loss[atom]) + ' ')

    log.debug('negative solution counter: ' + str(para.nega_count))
    log.debug('loss rejected counter: ' + str(para.loss_count))
    log.debug('delta rejected counter: ' + str(para.delta_count))
    log.info('------ Live long and prosper \\V/ ------')

`print_prog(var, para)`

Source code in src/vulcan/op.py

def print_prog(self, var, para):
    indx_max = np.nanargmax(para.where_varies_most)
    log.info(
        'Elapsed time: '
        + '{:.2e}'.format(var.t)
        + ' || Step number: '
        + str(para.count)
        + '/'
        + str(self.cfg.count_max)
    )
    log.info(
        'longdy = '
        + '{:.2e}'.format(var.longdy)
        + '      || longdy/dt = '
        + '{:.2e}'.format(var.longdydt)
        + '  || dt = '
        + '{:.2e}'.format(var.dt)
    )
    log.info('from nz = ' + str(int(indx_max / ni)) + ' and ' + species[indx_max % ni])
    log.info(spacer)

`save_out(var, atm, para)`

Source code in src/vulcan/op.py

def save_out(self, var, atm, para):
    output_dir, out_name = self.cfg.output_dir, self.cfg.out_name
    output_file = output_dir + out_name

    if not os.path.exists(output_dir):
        log.debug('The output directory assigned in vulcan_cfg.py does not exist.')
        log.debug(f'Directory {output_dir} created.')
        os.mkdir(output_dir)

    # convert lists into numpy arrays
    for key in var.var_evol_save:
        as_nparray = np.array(getattr(var, key))
        setattr(var, key, as_nparray)

    # plotting
    if self.cfg.use_plot_evo:
        self.plot_evo(var, atm)
    if self.cfg.use_plot_end:
        self.plot_update(
            var, atm, para, fpath=os.path.join(self.cfg.plot_dir, f'mix.{PLT_FMT}')
        )
    plt.close('all')

    # making the save dict
    var_save = {'species': species, 'nr': nr}

    for key in var.var_save:
        var_save[key] = getattr(var, key)
    if self.cfg.save_evolution:
        # slicing time-sequential data to reduce ouput filesize
        fq = self.cfg.save_evo_frq
        for key in var.var_evol_save:
            as_nparray = getattr(var, key)[::fq]
            setattr(var, key, as_nparray)
            var_save[key] = getattr(var, key)

    if self.cfg.output_humanread:
        # human-readable form, less efficient
        with open(output_file + '.var', 'w') as outfile:
            for var in var_save.keys():
                outfile.write(f'# {var} \n')
                outfile.write(str(var_save[var]) + '\n')

        with open(output_file + '.atm', 'w') as outfile:
            for var in vars(atm).keys():
                outfile.write(f'# {var} \n')
                outfile.write(str(vars(atm)[var]) + '\n')

        with open(output_file + '.par', 'w') as outfile:
            for var in vars(para).keys():
                outfile.write(f'# {var} \n')
                outfile.write(str(vars(para)[var]) + '\n')
    else:
        # pickled form, less safe
        from pickle import dump

        with open(output_file, 'wb') as outfile:
            dump(
                {'variable': var_save, 'atm': vars(atm), 'parameter': vars(para)},
                outfile,
                protocol=4,
            )

vulcan.op

op

ReadRate(vulcan_cfg)

cfg = vulcan_cfg instance-attribute

i = 1 instance-attribute

list_tri = [] instance-attribute

lim_lowT_rates(var, atm)

make_bins_read_cross(var, atm)

read_rate(var, atm)

remove_rate(var)

rev_rate(var, atm)

Integration(odesolver, output, vulcan_cfg)

all the operators required in the continuity equation: dn/dt + dphi/dz = P - L

or class incorporating the esential numerical operations?

atol = self.cfg.atol instance-attribute

cfg = vulcan_cfg instance-attribute

condense_sp_index = [(species.index(sp)) for sp in (self.cfg.condense_sp)] instance-attribute

mtol = self.cfg.mtol instance-attribute

non_gas_sp = self.cfg.non_gas_sp instance-attribute

non_gas_sp_index = [(species.index(sp)) for sp in (self.non_gas_sp)] instance-attribute

odesolver = odesolver instance-attribute

output = output instance-attribute

run_agni = run_agni instance-attribute

update_photo_frq = self.cfg.ini_update_photo_frq instance-attribute

use_settling = self.cfg.use_settling instance-attribute

backup(var)

conden(var, atm)

conv(var, para, atm, out=False, print_freq=100)

f_dy(var, para)

h2o_conden_evap_relax(var, atm)

h2o_conden_relax(var, atm)

nh3_conden_evap_relax(var, atm)

save_step(var, para)

stop(var, para, atm)

update_mu_dz(var, atm, make_atm)

update_phi_esc(var, atm)

ODESolver(vulcan_cfg)

atol = self.cfg.atol instance-attribute

cfg = vulcan_cfg instance-attribute

condense_sp_index = [(species.index(sp)) for sp in (self.cfg.condense_sp)] instance-attribute

fix_sp_bot_index = [(species.index(sp)) for sp in (self.cfg.use_fix_sp_bot.keys())] instance-attribute

fix_sp_bot_mix = np.array([(self.cfg.use_fix_sp_bot[sp]) for sp in (self.cfg.use_fix_sp_bot.keys())]) instance-attribute

mtol = self.cfg.mtol instance-attribute

non_gas_sp = self.cfg.non_gas_sp instance-attribute

non_gas_sp_index = [(species.index(sp)) for sp in (self.non_gas_sp)] instance-attribute

clip(var, para, atm, pos_cut=0, nega_cut=-1)

compute_J(var, atm)

compute_Jion(var, atm)

compute_flux(var, atm)

compute_tau(var, atm)

diffdf(y, atm)

diffdf_no_mol(y, atm)

diffdf_settling(y, atm)

diffdf_settling_vm(y, atm)

jac_tot(var, atm)

lhs_jac_fix_all_bot(var, atm)

lhs_jac_no_mol(var, atm)

lhs_jac_no_mol_fix_all_bot(var, atm)

lhs_jac_settling(var, atm)

lhs_jac_settling_vm(var, atm)

lhs_jac_tot(var, atm)

lhs_jac_tot_vm(var, atm)

loss(data_var)

print_lossBig(para)

print_nega(data_var, data_para)

reset_y(var, dt_reduc=0.5)

step_ok(var, para, loss_eps=0.1, rtol=0.6)

step_reject(var, para, loss_eps=0.1, rtol=0.6)

thomas_vec(a, b, c, d)

Ros2(vulcan_cfg)

naming_solver(para)

one_step(var, atm, para)

solver(var, atm, para)

solver_fix_all_bot(var, atm, para)

step_size(var, para, dt_var_min=2, dt_var_max=0.5, dt_min=1e-10, dt_max=1e+18)

store_bandM(a, nb, nn)

Output(vulcan_cfg)

cfg = vulcan_cfg instance-attribute

plot_TP(atm)

plot_evo(var, atm, plot_j=-1, plot_ymin=1e-16, dn=1)

`op`

`ReadRate(vulcan_cfg)`

`cfg = vulcan_cfg` `instance-attribute`

`i = 1` `instance-attribute`

`list_tri = []` `instance-attribute`

`lim_lowT_rates(var, atm)`

`make_bins_read_cross(var, atm)`

`read_rate(var, atm)`

`remove_rate(var)`

`rev_rate(var, atm)`

`Integration(odesolver, output, vulcan_cfg)`

`atol = self.cfg.atol` `instance-attribute`

`cfg = vulcan_cfg` `instance-attribute`

`condense_sp_index = [(species.index(sp)) for sp in (self.cfg.condense_sp)]` `instance-attribute`

`mtol = self.cfg.mtol` `instance-attribute`

`non_gas_sp = self.cfg.non_gas_sp` `instance-attribute`

`non_gas_sp_index = [(species.index(sp)) for sp in (self.non_gas_sp)]` `instance-attribute`

`odesolver = odesolver` `instance-attribute`

`output = output` `instance-attribute`

`run_agni = run_agni` `instance-attribute`

`update_photo_frq = self.cfg.ini_update_photo_frq` `instance-attribute`

`use_settling = self.cfg.use_settling` `instance-attribute`

`backup(var)`

`conden(var, atm)`

`conv(var, para, atm, out=False, print_freq=100)`

`f_dy(var, para)`

`h2o_conden_evap_relax(var, atm)`

`h2o_conden_relax(var, atm)`

`nh3_conden_evap_relax(var, atm)`

`save_step(var, para)`

`stop(var, para, atm)`

`update_mu_dz(var, atm, make_atm)`

`update_phi_esc(var, atm)`

`ODESolver(vulcan_cfg)`

`atol = self.cfg.atol` `instance-attribute`

`cfg = vulcan_cfg` `instance-attribute`

`condense_sp_index = [(species.index(sp)) for sp in (self.cfg.condense_sp)]` `instance-attribute`

`fix_sp_bot_index = [(species.index(sp)) for sp in (self.cfg.use_fix_sp_bot.keys())]` `instance-attribute`

`fix_sp_bot_mix = np.array([(self.cfg.use_fix_sp_bot[sp]) for sp in (self.cfg.use_fix_sp_bot.keys())])` `instance-attribute`

`mtol = self.cfg.mtol` `instance-attribute`

`non_gas_sp = self.cfg.non_gas_sp` `instance-attribute`

`non_gas_sp_index = [(species.index(sp)) for sp in (self.non_gas_sp)]` `instance-attribute`

`clip(var, para, atm, pos_cut=0, nega_cut=-1)`

`compute_J(var, atm)`

`compute_Jion(var, atm)`

`compute_flux(var, atm)`

`compute_tau(var, atm)`

`diffdf(y, atm)`

`diffdf_no_mol(y, atm)`

`diffdf_settling(y, atm)`

`diffdf_settling_vm(y, atm)`

`jac_tot(var, atm)`

`lhs_jac_fix_all_bot(var, atm)`

`lhs_jac_no_mol(var, atm)`

`lhs_jac_no_mol_fix_all_bot(var, atm)`

`lhs_jac_settling(var, atm)`

`lhs_jac_settling_vm(var, atm)`

`lhs_jac_tot(var, atm)`

`lhs_jac_tot_vm(var, atm)`

`loss(data_var)`

`print_lossBig(para)`

`print_nega(data_var, data_para)`

`reset_y(var, dt_reduc=0.5)`

`step_ok(var, para, loss_eps=0.1, rtol=0.6)`

`step_reject(var, para, loss_eps=0.1, rtol=0.6)`

`thomas_vec(a, b, c, d)`

`Ros2(vulcan_cfg)`

`naming_solver(para)`

`one_step(var, atm, para)`

`solver(var, atm, para)`

`solver_fix_all_bot(var, atm, para)`

`step_size(var, para, dt_var_min=2, dt_var_max=0.5, dt_min=1e-10, dt_max=1e+18)`

`store_bandM(a, nb, nn)`

`Output(vulcan_cfg)`

`cfg = vulcan_cfg` `instance-attribute`

`plot_TP(atm)`

`plot_evo(var, atm, plot_j=-1, plot_ymin=1e-16, dn=1)`

`plot_evo_inter(var, atm, plot_j=-1, plot_ymin=1e-16, dn=1)`

`plot_flux_update(var, atm, para)`

`plot_update(var, atm, para, fpath=None)`