6
\$\begingroup\$

Introduction

I'm new to python and made this code which computes the decay of sulphur dioxide (SO2) pollutant in a chemical reactor by ozone (O3) when it is hit by UV light.

Code description

Short explanation of the code: using the Cantera python package I define a set of properties for the gas (air + SO2 + O3) and the reactor. The reactions are stored in an external file (photolysis.cti) which I included after the code. Then the reactions are allowed to play out but at each time instant the rate of the reaction controlled by UV light is updated. This is done based on the instantaneous concentration of the species present and the UV light distribution. This can be modelled using several models I implemented. Some models (MPSS, MSSS) are more accurate but also more complex so the simulation will take longer if they are selected (~45s vs ~1s for the rest). After the simulation finished, a plot with the decay of SO2 with time is shown:

Text

Aim of this review

My aim is to improve both the python aspect (syntax, making it easier to read, more idiomatic, object oriented...) as well as the performance of the code. Also having the whole code in one single file feels not ideal. Any tips on how to structure the project are very welcome. I know the code is long but I tried to comment all functions to improve readability and made a very simple log to know what the code is doing.

Python code

The code I wrote is the following:

"""
Code to compute a CSTR UV reactor with Cantera
"""
import re
import time
import pandas as pd
import numpy as np
import cantera as ct
import discretize.CylMesh
import matplotlib.pyplot as plt
from scipy.constants import N_A, h, c


def reactor_volume(L, r_outer, r_inner):
    """
    Computes the volume of an annular reactor (hollow cylinder).

    Parameters [units]
    ------------------
    L : int, float
        Length of the reactor/lamp [m]
    r_outer : int, float
        Outer radius [m]
    r_outer : int, float
        Inner radius [m]

    Returns [units]
    ---------------
    V_r : float
        The reactor volume [m**2]
    """
    return np.pi*L*(r_outer**2-r_inner**2)


def make_cyl_mesh(L, r, ncr, ncz):
    """
    Makes a cylindrical mesh object with one cell in the azimuthal direction.

    Parameters [units]
    ------------------
    L : int, float
        Length of the reactor/lamp [m]
    r : int, float
        Cylinder radius [m]
    ncr : int
        Number of cells in the radial direction
    ncz : int
        Number of cells in the axial direction

    Returns [units]
    ---------------
    mesh : obj
        The mesh cylinder object
    """
    dr = r/ncr              # cell width r
    dz = L/ncz              # cell width z
    hr = dr*np.ones(ncr)
    hz = dz*np.ones(ncz)
    offset = [0, 0, -L/2]   # offset to center z on half the lamp length
    return discretize.CylMesh([hr, 2*np.pi, hz], offset)


def sigma_to_eps(sigma):
    """
    Converts the molar absorbtion cross section of a gas to the molar
    absorbtion coefficient in S.I. units.

    Parameters [units]
    ------------------
    sigma : int, float
        Absorbtion cross section (log base e) [cm**2]

    Returns [units]
    ---------------
    eps : float
        molar absorbtion coefficient (log base 10) [m**2 mol**-1]
    """
    return (sigma*N_A)/(1e4*np.log(10))


def ppm_to_C(ppm):
    """
    Converts ppm (parts per million) concentration of a gas to
    concentration in S.I. units.

    Parameters [units]
    ------------------
    ppm : int, float
        ppm concentration of a gas [-]

    Returns [units]
    ---------------
    C : float
        concentration [mol m**-3]
    """
    rho_air = 1.2  # [kg m**-3]
    Mr_air = 28.96  # [g mol**-1] Note: this is for DRY air
    Vm_air = molar_volume(rho_air, Mr_air)  # [m**3 mol**-1]
    C = ppm/(Vm_air*1e6)
    return C


def molar_volume(rho, Mr):
    """
    Computes the molar volume of a gas.

    Parameters [units]
    ------------------
    rho : int, float
        density [kg m**-3]
    Mr : int, float
        molar mass [g mol**-1]

    Returns [units]
    ---------------
    Vm : float
        molar volume [m**3 mol**-1]
    """
    return (Mr/rho)*1e-3


def exp_medium_absorption(r, ppm=300):
    """
    Computes the absorption of the medium considered in the fluence 
    rate models. Initially, I will treat only SO2 as absorbing species.

    Parameters [units]
    ------------------
    r : int, float, array
        Radial position (i.e. of a mesh element) [m]
    ppm : int, float
        ppm concentration of a gas [-]

    Returns [units]
    ---------------
    ema : int, float, array
        Exponential medium absorption [-]
    """
    sigma = 1.551258e-19  # Value for SO2 [cm**2]
    sigma_air = sigma
    eps_air = sigma_to_eps(sigma_air)
    C_air = ppm_to_C(ppm)
    return np.exp(-np.log(10)*eps_air*C_air*r)


def coord_ls(L, nls):
    """
    Computes the coordinates of the number light source elements 
    which have been chosen to model the UV lamp.

    Parameters [units]
    ------------------
    L : int, float
        Length of the reactor/lamp [m]
    nls : int
        Number of light source elements

    Returns [units]
    ---------------
    coord_ls : array
        Coordinates of the number light source elements [m]
    """
    d_ls = L/nls
    return np.arange(-L/2, L/2+d_ls, d_ls)


def optimize_nls(case_variables, tol=1E-2):
    """
    Computes the optimum number of light sources (nls) for which to
    model a UV lamp (for the MPSS or MSSS models). This optimum is 
    the lowest nls that makes the fluence rate computed by the model
    (without medium absorbtion) match the fluence rate computed by 
    the LSI model for a given tolerance. Based on Powell and Lawryshyn 
    2015.
    Note: this routine takes some time, proceed with caution.

    Parameters [units]
    ------------------
    case_variables : array
        Array contining the following case variables:
        case_variables = [P, L, eff, r_outer, r_inner, ncr, ncz]    
        P       = lamp power [W]
        L       = lamp length [m]
        eff     = lamp efficiency at 254 nm wavelength
        r_outer = reactor radius [m]
        r_inner = quartz sleeve radius [m]
        ncr     = number of mesh cells in r
        ncz     = number of mesh cells in z
    tol : int, float
        Tolerance for the optimization routine

    Returns [units]
    ---------------
    nls : int
        (Optimum) number of light source elements
    """
    P, eff, L, r_outer, r_inner, ncr, ncz = case_variables
    V_r = reactor_volume(L, r_outer, r_inner)
    mesh_outer = make_cyl_mesh(L, r_outer, ncr, ncz)
    mesh_inner = make_cyl_mesh(L, r_inner, ncr, ncz)
    err, nls = 100, 10  # Initialize error and smallest possible number of light sources
    # Compute reference LSI vafr
    I_LSI = compute_vafr(case_variables, model='LSI')
    # Loop to find the optimum value for the number of light sources
    while err > tol*I_LSI:
        nls += 1
        I_i_outer = per_cell_fluence_rate_MPSS(
            mesh_outer, P, eff, L, nls=nls, no_abs=True)
        I_i_inner = per_cell_fluence_rate_MPSS(
            mesh_inner, P, eff, L, nls=nls, no_abs=True)
        I_MPSS = volume_avg_fluence_rate(
            I_i_outer, mesh_outer, I_i_inner, mesh_inner, V_r)
        err = np.abs(I_MPSS-I_LSI)
    return nls


def per_cell_fluence_rate_RAD(mesh, P, eff, L, no_abs=False, ppm=300):
    """
    Computes the fluence rate for each cell in a given mesh following
    the RAD model (Coenen 2013).

    Parameters [units]
    ------------------
    mesh : obj
        Mesh object
    P : int, float
        Lamp power [W]
    eff : float
        Lamp efficiency at 254 nm wavelength
    L : int, float
        Length of the reactor/lamp [m]
    no_abs : bool
        Flag to enable/disable medium absorption
    ppm : int, float
        ppm concentration of a gas [-]

    Returns [units]
    ---------------
    I_i : array
        Fluence rate for each cell in the mesh [W m**-2]
    """
    R = mesh.gridCC[:, 0]
    if no_abs == True:
        I_i = ((P*eff)/(2*np.pi*R*L))
    else:
        I_i = ((P*eff)/(2*np.pi*R*L))*exp_medium_absorption(R, ppm=ppm)
    return I_i


def per_cell_fluence_rate_LSI(mesh, P, eff, L):
    """
    Computes the fluence rate for each cell in a given mesh following
    the LSI model (Blatchley 1998). It assumes no medium absorption

    Parameters [units]
    ------------------
    mesh : obj
        Mesh object
    P : int, float
        Lamp power [W]
    eff : float
        Lamp efficiency at 254 nm wavelength
    L : int, float
        Length of the reactor/lamp [m]

    Returns [units]
    ---------------
    I_i : array
        Fluence rate for each cell in the mesh [W m**-2]
    """
    R = mesh.gridCC[:, 0]
    H = mesh.gridCC[:, 2]
    I_i = ((P*eff)/(4*np.pi*L*R))*(np.arctan((L/2+H)/R)+np.arctan((L/2-H)/R))
    return I_i


def per_cell_fluence_rate_RAD_LSI(mesh, P, eff, L):
    """
    Computes the fluence rate for each cell in a given mesh following
    the RAD-LSI model (Liu 2004). It assumes no medium absorption

    Parameters [units]
    ------------------
    mesh : obj
        Mesh object
    P : int, float
        Lamp power [W]
    eff : float
        Lamp efficiency at 254 nm wavelength
    L : int, float
        Length of the reactor/lamp [m]

    Returns [units]
    ---------------
    I_i : array
        Fluence rate for each cell in the mesh [W m**-2]
    """
    R = mesh.gridCC[:, 0]
    H = mesh.gridCC[:, 2]
    I_i_LSI = ((P*eff)/(4*np.pi*L*R)) * \
        (np.arctan((L/2+H)/R)+np.arctan((L/2-H)/R))
    I_i_RAD = ((P*eff)/(2*np.pi*R*L))
    I_i = np.minimum(I_i_RAD, I_i_LSI)
    return I_i


def per_cell_fluence_rate_MPSS(mesh, P, eff, L, nls, no_abs=False, ppm=300):
    """
    Computes the fluence rate for each cell in a given mesh following
    the MPSS model (Bolton 2000).

    Parameters [units]
    ------------------
    mesh : obj
        Mesh object
    P : int, float
        Lamp power [W]
    eff : float
        Lamp efficiency at 254 nm wavelength
    L : int, float
        Length of the reactor/lamp [m]
    nls : int
        Number of light source elements
    no_abs : bool
        Flag to enable/disable medium absorption
    ppm : int, float
        ppm concentration of a gas [-]

    Returns [units]
    ---------------
    I_i : array
        Fluence rate for each cell in the mesh [W m**-2]
    """
    x_ls = coord_ls(L, nls)
    R = mesh.gridCC[:, 0]
    H = mesh.gridCC[:, 2]
    I_i = np.zeros(len(R))
    for i in range(len(R)):
        r = R[i]*np.ones(len(x_ls))
        h = H[i]*np.ones(len(x_ls))
        rho = np.sqrt(r**2+(h-x_ls)**2)
        if no_abs == True:
            I_ls = (((P*eff)/nls)/(4*np.pi*rho**2))
        else:
            I_ls = (((P*eff)/nls)/(4*np.pi*rho**2)) * \
                exp_medium_absorption(r, ppm)
        I_i[i] = sum(I_ls)
    return I_i


def per_cell_fluence_rate_MSSS(mesh, P, eff, L, nls, no_abs=False, ppm=300):
    """
    Computes the fluence rate for each cell in a given mesh following
    the MSSS model (Liu 2004).

    Parameters [units]
    ------------------
    mesh : obj
        Mesh object
    P : int, float
        Lamp power [W]
    eff : float
        Lamp efficiency at 254 nm wavelength
    L : int, float
        Length of the reactor/lamp [m]
    nls : int
        Number of light source elements
    no_abs : bool
        Flag to enable/disable medium absorption
    ppm : int, float
        ppm concentration of a gas [-]

    Returns [units]
    ---------------
    I_i : array
        Fluence rate for each cell in the mesh [W m**-2]
    """
    x_ls = coord_ls(L, nls)
    R = mesh.gridCC[:, 0]
    H = mesh.gridCC[:, 2]
    I_i = np.zeros(len(R))
    for i in range(len(R)):
        r = R[i]*np.ones(len(x_ls))
        h = H[i]*np.ones(len(x_ls))
        rho = np.sqrt(r**2+(h-x_ls)**2)
        if no_abs == True:
            I_ls = (((P*eff)/nls)/(4*np.pi*rho**2))
        else:
            I_ls = (((P*eff)/nls)/(4*np.pi*rho**2)) * \
                exp_medium_absorption(r, ppm)
        I_ls = I_ls*np.cos(np.arctan(r/h))
        I_i[i] = sum(I_ls)
    return I_i


def per_cell_fluence_rate(mesh, P, eff, L, model='RAD', nls=40, ppm=300):
    """
    Generic function that computes the fluence rate for each cell 
    in a given mesh for various model.

    Parameters [units]
    ------------------
    mesh : obj
        Mesh object
    P : int, float
        Lamp power [W]
    eff : float
        Lamp efficiency at 254 nm wavelength
    L : int, float
        Length of the reactor/lamp [m]
    model : str
        name of the model ('RAD', 'LSI', 'MPSS', 'MSSS' or 'RAD_LSI')
    nls : int
        Number of light source elements
    ppm : int, float
        ppm concentration of a gas [-]

    Returns [units]
    ---------------
    I_i : array
        Fluence rate for each cell in the mesh [W m**-2]
    """
    if isinstance(model, str) == 0:
        raise TypeError("The variable 'model' has to be a string")

    if model == 'RAD':
        I_i = per_cell_fluence_rate_RAD(mesh, P, eff, L, ppm=ppm)
    elif model == 'LSI':
        I_i = per_cell_fluence_rate_LSI(mesh, P, eff, L)
    elif model == 'MPSS':
        I_i = per_cell_fluence_rate_MPSS(mesh, P, eff, L, nls=nls, ppm=ppm)
    elif model == 'MSSS':
        I_i = per_cell_fluence_rate_MSSS(mesh, P, eff, L, nls=nls, ppm=ppm)
    elif model == 'RAD_LSI':
        I_i = per_cell_fluence_rate_RAD_LSI(mesh, P, eff, L)
    else:
        raise ValueError(
            "Only 'RAD', 'LSI', 'MPSS', 'MSSS' or 'RAD_LSI' models are implemented")
    return I_i


def volume_avg_fluence_rate(I_i_outer, mesh_outer, I_i_inner, mesh_inner, V_r):
    """
    Computes the volume averaged fluence rate (vafr) for the
    annular reactor (hollow cylinder) by subtracting the
    average of the small cylinder (_inner) from the large 
    cylinder (_outer).

    Parameters [units]
    ------------------
    I_i_outer : array
        Fluence rate for each cell in the outer cylinder mesh [W m**-2]
    mesh_outer : obj
        outer cylinder mesh object
    I_i_inner : array
        Fluence rate for each cell in the inner cylinder mesh [W m**-2]
    mesh_inner : obj
        Inner cylinder mesh object
    V_r : float
        The reactor volume [m**2]

    Returns [units]
    ---------------
    I_vol_avg : array
        Volume average fluence rate for the annular reactor [W m**-2]
    """
    # volume averaged fluence rate (vafr)
    I_outer = np.sum(I_i_outer*mesh_outer.vol)
    I_inner = np.sum(I_i_inner*mesh_inner.vol)
    return (I_outer - I_inner)/V_r


def compute_vafr(case_variables, model='RAD', nls=40, optimize=False, ppm=300):
    """
    Generic function to compute the volume average fluence rate.

    Parameters [units]
    ------------------
    case_variables : array
        Array contining the following case variables:
        case_variables = [P, L, eff, r_outer, r_inner, ncr, ncz]    
        P       = lamp power [W]
        L       = lamp length [m]
        eff     = lamp efficiency at 254 nm wavelength
        r_outer = reactor radius [m]
        r_inner = quartz sleeve radius [m]
        ncr     = number of mesh cells in r
        ncz     = number of mesh cells in z
    model : str
        name of the model ('RAD', 'LSI', 'MPSS', 'MSSS' or 'RAD_LSI')
    nls : int
        Number of light source elements
    no_abs : bool
        Flag to enable/disable nls optimization for the MPSS or MSSS models
    ppm : int, float
        ppm concentration of a gas [-]

    Returns [units]
    ---------------
    I_vol_avg : array
        Volume average fluence rate for the annular reactor [W m**-2]
    """
    # Unpack the case variables
    P, eff, L, r_outer, r_inner, ncr, ncz = case_variables

    # Compute reactor volume
    V_r = reactor_volume(L, r_outer, r_inner)

    # Make the outer and inner mesh objects
    mesh_outer = make_cyl_mesh(L, r_outer, ncr, ncz)
    mesh_inner = make_cyl_mesh(L, r_inner, ncr, ncz)

    # Find optimum nls if required
    if (model == 'MPSS' or model == 'MSSS') and optimize == True:
        nls = optimize_nls(case_variables, tol=1E-2)

    # Calculate the per-cell fluence rate for the inner and outer meshes
    I_i_outer = per_cell_fluence_rate(
        mesh_outer, P, eff, L, model=model, nls=nls, ppm=ppm)
    I_i_inner = per_cell_fluence_rate(
        mesh_inner, P, eff, L, model=model, nls=nls, ppm=ppm)

    # Compute volume average LSI fluence rate
    I_vol_avg = volume_avg_fluence_rate(
        I_i_outer, mesh_outer, I_i_inner, mesh_inner, V_r)
    return I_vol_avg


def set_gas(mechanism, T, P, C):
    """
    Define the properties of the cantera gas object.

    Parameters [units]
    ------------------
    mechanism : string
        Path to the CTI mechanism to be used
    T : int, float
        Temperature [K]
    P : int, float
        Pressure [in atm]
    C : dict
        Dictionary containing the concentrations of the different gas species [molar fraction]

    Returns [units]
    ---------------
    gas : obj
        cantera gas object
    """
    gas = ct.Solution('data/photolysis.cti')
    gas.TPX = T, P, C
    return gas


def set_reactor(gas, V_r, residence_t, p_valve_coeff):
    """
    Set-up the reactor network using Cantera functions and syntax.

    Parameters [units]
    ------------------
    gas : obj
        Cantera gas object
    V_r : float
        The reactor volume [m**2]
    residence_t : int, float
        Time that the gas will stay in the reactor [s]
    p_valve_coeff : int, float
        This is the "conductance" of the pressure valve and will 
        determine its efficiency in holding the reactor pressure 
        to the desired conditions. Set to 0.01 for default

    Returns [units]
    ---------------
    stirred_reactor : obj
        Cantera reactor object
    reactor_network : obj
        Cantera network object
    """
    fuel_air_mixture_tank = ct.Reservoir(gas)
    exhaust = ct.Reservoir(gas)
    stirred_reactor = ct.IdealGasReactor(gas, energy='off', volume=V_r)
    mass_flow_controller = ct.MassFlowController(
        upstream=fuel_air_mixture_tank, downstream=stirred_reactor, mdot=stirred_reactor.mass/residence_t)
    pressure_regulator = ct.Valve(
        upstream=stirred_reactor, downstream=exhaust, K=p_valve_coeff)
    reactor_network = ct.ReactorNet([stirred_reactor])
    return [stirred_reactor, reactor_network]


def intialize_results(stirred_reactor):
    """
    Initializes the DataFrame where results will be stored.

    Parameters [units]
    ------------------
    stirred_reactor : obj
        Cantera reactor object

    Returns [units]
    ---------------
    time_history : DataFrame
        Pandas DataFrame to store time evolution of reasults. It
        includes the gas's mass, volume, temperature and species
        concentration
    """
    column_names = [stirred_reactor.component_name(
        item) for item in range(stirred_reactor.n_vars)]
    time_history = pd.DataFrame(columns=column_names)
    return time_history


def update_results(time_history, stirred_reactor, t):
    """
    Updates the results at a given time instant.

    Parameters [units]
    ------------------
    time_history : DataFrame
        Pandas DataFrame to store time evolution of reasults. It
        includes the gas's mass, volume, temperature and species
        concentration
    stirred_reactor : obj
        Cantera reactor object
    t : int, float
        Current simulation time [s]

    Returns [units]
    ---------------
    time_history : DataFrame
        Pandas DataFrame to store time evolution of reasults. It
        includes the gas's mass, volume, temperature and species
        concentration
    """
    state = np.hstack([stirred_reactor.mass, stirred_reactor.volume,
                       stirred_reactor.T, stirred_reactor.thermo.X])
    time_history.loc[t] = state
    return time_history


def uv_update_reactions(stirred_reactor, case_variables, model='RAD', nls=40):
    """
    Modifies the rate constant of the photolysis reaction of ozone. This
    has to be done for each loop iteration as the concentrations of
    the different gas species are changing, which will affect the absorption
    of the medium.

    Parameters [units]
    ------------------
    stirred_reactor : obj
        Cantera reactor object
    case_variables : array
        Array contining the following case variables:
        case_variables = [P, L, eff, r_outer, r_inner, ncr, ncz]    
        P       = lamp power [W]
        L       = lamp length [m]
        eff     = lamp efficiency at 254 nm wavelength
        r_outer = reactor radius [m]
        r_inner = quartz sleeve radius [m]
        ncr     = number of mesh cells in r
        ncz     = number of mesh cells in z
    model : str
        name of the model ('RAD', 'LSI', 'MPSS', 'MSSS' or 'RAD_LSI')
    nls : int
        Number of light source elements

    Returns [units]
    ---------------
    None
    """
    ppm = stirred_reactor.thermo.X[stirred_reactor.component_index(
        'SO2')-3]*1e6
    k_uv = (0.9) * ((254E-9)/(N_A*h*c)) * (np.log(10)) * sigma_to_eps(1.132935E-17) * \
        compute_vafr(case_variables, model=model, nls=nls, ppm=ppm)
    ID = np.size(gas.reactions())-1
    reaction = gas.reactions()[ID]
    reaction.rate = ct.Arrhenius(A=k_uv, b=0, E=0)
    gas.modify_reaction(ID, reaction)
    return None


def time_loop(reactor_network, stirred_reactor, case_variables, time_history, max_simulation_t, model='RAD'):
    """
    Computes the time loop iterations for the chemical reactor
    until the maximum simulation time is reached.

    Parameters [units]
    ------------------
    stirred_reactor : obj
        Cantera reactor object
    reactor_network : obj
        Cantera network object
    case_variables : array
        Array contining the following case variables:
        case_variables = [P, L, eff, r_outer, r_inner, ncr, ncz]    
        P       = lamp power [W]
        L       = lamp length [m]
        eff     = lamp efficiency at 254 nm wavelength
        r_outer = reactor radius [m]
        r_inner = quartz sleeve radius [m]
        ncr     = number of mesh cells in r
        ncz     = number of mesh cells in z
    time_history : DataFrame
        Pandas DataFrame to store time evolution of reasults. It
        includes the gas's mass, volume, temperature and species
        concentration
    max_simulation_t : int
        Maximum simulation time [s]
    model : str
        name of the model ('RAD', 'LSI', 'MPSS', 'MSSS' or 'RAD_LSI')

    Returns [units]
    ---------------
    time_history : DataFrame
        Pandas DataFrame to store time evolution of reasults. It
        includes the gas's mass, volume, temperature and species
        concentration
    """
    counter, t = 1, 0
    while t < max_simulation_t:
        if (t == 0) or (counter % 10 == 0): # Update every 10 iterations
            if (model == 'MPSS') or (model == 'MSSS'):
                if t == 0:
                    nls = optimize_nls(case_variables, tol=1E-2)
                uv_update_reactions(
                    stirred_reactor, case_variables, model=model, nls=nls)
            else:
                uv_update_reactions(stirred_reactor, case_variables, model=model)
            time_history = update_results(time_history, stirred_reactor, t)
        t = reactor_network.step()
        counter += 1
    return time_history


def run_simulation(gas, V_r, residence_t, p_valve_coeff, case_variables, max_simulation_t, model='RAD'):
    """
    Generic function to run the chemical reactor with UV light
    modeling.

    Parameters [units]
    ------------------
    gas : obj
        Cantera gas object
    V_r : float
        The reactor volume [m**2]
    residence_t : int, float
        Time that the gas will stay in the reactor [s]
    p_valve_coeff : int, float
        This is the "conductance" of the pressure valve and will 
        determine its efficiency in holding the reactor pressure 
        to the desired conditions. Set to 0.01 for default
    case_variables : array
        Array contining the following case variables:
        case_variables = [P, L, eff, r_outer, r_inner, ncr, ncz]    
        P       = lamp power [W]
        L       = lamp length [m]
        eff     = lamp efficiency at 254 nm wavelength
        r_outer = reactor radius [m]
        r_inner = quartz sleeve radius [m]
        ncr     = number of mesh cells in r
        ncz     = number of mesh cells in z
    max_simulation_t : int
        Maximum simulation time [s]
    model : str
        name of the model ('RAD', 'LSI', 'MPSS', 'MSSS' or 'RAD_LSI')

    Returns [units]
    ---------------
    time_history : DataFrame
        Pandas DataFrame to store time evolution of reasults. It
        includes the gas's mass, volume, temperature and species
        concentration
    """
    # Define the reactor object
    [stirred_reactor, reactor_network] = set_reactor(
        gas, V_r, residence_t, p_valve_coeff)
    # Initialize the results dataFrame
    time_history = intialize_results(stirred_reactor)
    # Run the time loop
    time_history = time_loop(reactor_network, stirred_reactor, case_variables, time_history,
                             max_simulation_t,  model=model)
    return time_history


def plot_results(specie, log=False):
    """
    Plots the results from the reactor simulation.

    Parameters [units]
    ------------------
    specie : str
        Specie for which to plot the time evolution
    log : bool
        Flag to enable/disable the x axis as logarithmic

    Returns [units]
    ---------------
    None
    """
    plt.style.use('ggplot')
    plt.style.use('seaborn-pastel')
    plt.rcParams['axes.labelsize'] = 18
    plt.rcParams['xtick.labelsize'] = 14
    plt.rcParams['ytick.labelsize'] = 14
    plt.rcParams['figure.autolayout'] = True

    plt.figure()
    if log == True:
        plt.semilogx(time_history.index, time_history[specie]*1e6, '-')
    else:
        plt.plot(time_history.index, time_history[specie]*1e6, '-')
    plt.xlabel('Time (s)')
    plt.ylabel(re.sub(r"([0-9]+(\.[0-9]+)?)", r"_\1",
                      '$'+specie+'$' + r' Mole Fraction : $ppmv$'))
    plt.show()
    return None


if __name__ == "__main__":
    start = time.time()

    print('Defining the gas properties')
    gas = set_gas('photolysis.cti', 295, ct.one_atm,
                  {'N2': 0.78014, 'O2': 0.20946, 'O3': 0.0001, 'SO2': 0.0003, 'H2O': 0.01})

    print('Defining the geometry variables')
    case_variables = P, eff, L, r_outer, r_inner, ncr, ncz = [
        17, 0.33, 0.28, 0.0425, 0.0115, 20, 20]
    V_r = reactor_volume(L, r_outer, r_inner)

    print('Defining the simulation variables')
    p_valve_coeff, max_p_rise, residence_t, max_simulation_t = 0.01, 0.01, 10, 100

    print('Running the simulation')
    time_history = run_simulation(gas, V_r, residence_t, p_valve_coeff, case_variables,
                                  max_simulation_t, model='RAD')

    end = time.time()

    print('Plotting the simulation results')
    plot_results('SO2', log=True)

    print(f'Elapsed time: {end - start:.3f}s')

Kinetic mechanism

If you wish to run the code you will require the kinetic mechanism. Save it as photolysis.cti:

"""
-----------------------------------------------------------------------|
                                                                       |
                         Photolysis mechanism                          |
                                                                       |
-----------------------------------------------------------------------|
"""

units(length='cm', time='s', quantity='molec', act_energy='K')

ideal_gas(name='gas',
          elements="O H N S",
          species="""O      O(1D)  O2     O3     H2O    HO2    OH     N2
                     SO2    SO3    HSO3   H2SO4  H""",
          reactions='all',
          initial_state=state(temperature=300.0, pressure=OneAtm))

#-------------------------------------------------------------------------------
# Species data
#-------------------------------------------------------------------------------

species(name='O',
        atoms='O:1',
        thermo=(NASA([200.00, 1000.00],
                     [ 2.94642900E+00, -1.63816600E-03,  2.42103200E-06,
                      -1.60284300E-09,  3.89069600E-13,  2.91476400E+04,
                       2.96399500E+00]),
                NASA([1000.00, 5000.00],
                     [ 2.54206000E+00, -2.75506200E-05, -3.10280300E-09,
                       4.55106700E-12, -4.36805200E-16,  2.92308000E+04,
                       4.92030800E+00])),
        note='120186')

species(name='O(1D)',
        atoms='O:1',
        thermo=(NASA([200.00, 1000.00],
                     [ 2.49993786E+00,  1.71935346E-07, -3.45215267E-10,
                       3.71342028E-13, -1.70964494E-16,  5.19965317E+04,
                       4.61684555E+00]),
                NASA([1000.00, 6000.00],
                     [ 2.49368475E+00,  1.37617903E-05, -1.00401058E-08,
                       2.76012182E-12, -2.01597513E-16,  5.19986304E+04,
                       4.65050950E+00])),
        note='121286')

species(name='O2',
        atoms='O:2',
        thermo=(NASA([200.00, 1000.00],
                     [ 3.21293600E+00,  1.12748600E-03, -5.75615000E-07,
                       1.31387700E-09, -8.76855400E-13, -1.00524900E+03,
                       6.03473800E+00]),
                NASA([1000.00, 5000.00],
                     [ 3.69757800E+00,  6.13519700E-04, -1.25884200E-07,
                       1.77528100E-11, -1.13643500E-15, -1.23393000E+03,
                       3.18916600E+00])),
        note='121386')

species(name='O3',
        atoms='O:3',
        thermo=(NASA([200.00, 1000.00],
                     [ 2.46260900E+00,  9.58278100E-03, -7.08735900E-06,
                       1.36336800E-09,  2.96964700E-13,  1.60615200E+04,
                       1.21418700E+01]),
                NASA([1000.00, 5000.00],
                     [ 5.42937100E+00,  1.82038000E-03, -7.70560700E-07,
                       1.49929300E-10, -1.07556300E-14,  1.52352700E+04,
                      -3.26638700E+00])),
        note='121286')

species(name='H2O',
        atoms='H:2 O:1',
        thermo=(NASA([200.00, 1000.00],
                     [ 3.38684200E+00,  3.47498200E-03, -6.35469600E-06,
                       6.96858100E-09, -2.50658800E-12, -3.02081100E+04,
                       2.59023300E+00]),
                NASA([1000.00, 5000.00],
                     [ 2.67214600E+00,  3.05629300E-03, -8.73026000E-07,
                       1.20099600E-10, -6.39161800E-15, -2.98992100E+04,
                       6.86281700E+00])),
        note='20387')

species(name='HO2',
        atoms='H:1 O:2',
        thermo=(NASA([200.00, 1000.00],
                     [ 2.97996300E+00,  4.99669700E-03, -3.79099700E-06,
                       2.35419200E-09, -8.08902400E-13,  1.76227400E+02,
                       9.22272400E+00]),
                NASA([1000.00, 5000.00],
                     [ 4.07219100E+00,  2.13129600E-03, -5.30814500E-07,
                       6.11226900E-11, -2.84116500E-15, -1.57972700E+02,
                       3.47602900E+00])),
        note='20387')

species(name='OH',
        atoms='H:1 O:1',
        thermo=(NASA([200.00, 1000.00],
                     [ 3.63726600E+00,  1.85091000E-04, -1.67616500E-06,
                       2.38720300E-09, -8.43144200E-13,  3.60678200E+03,
                       1.35886000E+00]),
                NASA([1000.00, 5000.00],
                     [ 2.88273000E+00,  1.01397400E-03, -2.27687700E-07,
                       2.17468400E-11, -5.12630500E-16,  3.88688800E+03,
                       5.59571200E+00])),
        note='121286')

species(name='N2',
        atoms='N:2',
        thermo=(NASA([100.00, 1000.00],
                     [ 3.29867700E+00,  1.40824000E-03, -3.96322200E-06,
                       5.64151500E-09, -2.44485500E-12, -1.02090000E+03,
                       3.95037200E+00]),
                NASA([1000.00, 5000.00],
                     [ 2.92664000E+00,  1.48797700E-03, -5.68476100E-07,
                       1.00970400E-10, -6.75335100E-15, -9.22797700E+02,
                       5.98052800E+00])),
        note='121286')

species(name='SO2',
        atoms='O:2 S:1',
        thermo=(NASA([200.00, 1000.00],
                     [ 2.91143900E+00,  8.10302200E-03, -6.90671000E-06,
                       3.32901600E-09, -8.77712100E-13, -3.68788200E+04,
                       1.11174000E+01]),
                NASA([1000.00, 5000.00],
                     [ 5.25449800E+00,  1.97854500E-03, -8.20422600E-07,
                       1.57638300E-10, -1.12045100E-14, -3.75688600E+04,
                      -1.14605600E+00])),
        note='121286')

species(name='SO3',
        atoms='O:3 S:1',
        thermo=(NASA([200.00, 1000.00],
                     [ 2.57528300E+00,  1.51509200E-02, -1.22987200E-05,
                       4.24025700E-09, -5.26681200E-13, -4.89441100E+04,
                       1.21951200E+01]),
                NASA([1000.00, 5000.00],
                     [ 7.05066800E+00,  3.24656000E-03, -1.40889700E-06,
                       2.72153500E-10, -1.94236500E-14, -5.02066800E+04,
                      -1.10644300E+01])),
        note='121286')

species(name='HSO3',
        atoms='H:1 O:3 S:1',
        thermo=(NASA([200.00, 1000.00],
                     [ 3.18471474E+00,  1.92277645E-02, -1.85975873E-05,
                       7.95502482E-09, -9.42254910E-13, -4.42156599E+04,
                       1.30496771E+01]),
                NASA([1000.00, 6000.00],
                     [ 8.19420431E+00,  3.77828016E-03, -1.34903219E-06,
                       2.17197023E-10, -1.29874848E-14, -4.55013223E+04,
                      -1.23824851E+01])),
        note='HO-SO2T10/10')

species(name='H2SO4',
        atoms='H:2 O:4 S:1',
        thermo=(NASA([200.00, 1000.00],
                     [ 4.53388173E+00,  3.10347679E-02, -4.10421795E-05,
                       2.95752341E-08, -8.81459071E-12, -9.05459072E+04,
                       3.93961412E+00]),
                NASA([1000.00, 6000.00],
                     [ 1.13355392E+01,  5.60829109E-03, -1.94574192E-06,
                       3.07136054E-10, -1.81109544E-14, -9.21087435E+04,
                      -2.96094003E+01])),
        note='T8/03')

species(name='H',
        atoms='H:1',
        thermo=(NASA([200.00, 1000.00],
                     [ 2.50000000E+00,  0.00000000E+00,  0.00000000E+00,
                       0.00000000E+00,  0.00000000E+00,  2.54716300E+04,
                      -4.60117600E-01]),
                NASA([1000.00, 5000.00],
                     [ 2.50000000E+00,  0.00000000E+00,  0.00000000E+00,
                       0.00000000E+00,  0.00000000E+00,  2.54716300E+04,
                      -4.60117600E-01])),
        note='120186')

#-------------------------------------------------------------------------------
# Reaction data
#-------------------------------------------------------------------------------

# Reaction 1
reaction('O + O3 => 2 O2', [8.000000e+12, 0.0, 2060.0])

# Reaction 2
reaction('O(1D) + H2O => 2 OH', [1.630000e+10, 0.0, 65.0])

# Reaction 3
reaction('O + OH => O2 + H', [2.400000e+11, 0.0, 110.0])

# Reaction 4
reaction('O + HO2 => OH + O2', [2.700000e+11, 0.0, 224.0])

# Reaction 5
three_body_reaction('H + O2 + M => HO2 + M', [2.170000e+29, 1.1, 0.0])

# Reaction 6
reaction('OH + O3 => HO2 + O2', [1.700000e+12, 0.0, 940.0])

# Reaction 7
reaction('OH + HO2 => H2O + O2', [4.800000e+11, 0.0, 250.0])

# Reaction 8
reaction('OH + SO2 => HSO3', [7.070000e+11, 0.7, 0.0])

# Reaction 9
three_body_reaction('O + SO2 + M => SO3 + M', [4.000000e+32, 0.0, 1000.0])

# Reaction 10
reaction('SO3 + H2O => H2SO4', [1.200000e+15, 0.0, 0.0])

# Reaction 11
reaction('HSO3 + O2 => HO2 + SO3', [1.300000e+12, 0.0, 330.0])

#-------------------------------------------------------------------------------
# Photolysis reaction
#-------------------------------------------------------------------------------

# Reaction 12
reaction('O3 => O(1D) + O2', [1.000000e+00, 0.0, 0.0])
\$\endgroup\$

1 Answer 1

4
\$\begingroup\$

Overall

This code is fairly well-structured and very well-documented. Excellent job.

Type hints

It's good that you specified types here:

    L : int, float
        Length of the reactor/lamp [m]
    r_outer : int, float
        Outer radius [m]
    r_outer : int, float

Better would be to codify these in type hints:

def reactor_volume(
    L: Union[int, float], 
    r_outer: Union[int, float], 
    r_inner: Union[int, float],
) -> float:

Some alternatives to the union are explained here:

https://stackoverflow.com/questions/50928592/mypy-type-hint-unionfloat-int-is-there-a-number-type

Order of operations

(sigma*N_A)/(1e4*np.log(10))

can be

sigma * N_A * 1e-4 / np.log(10)

Simplification

The exp will cancel out the log here:

np.exp(-np.log(10)*eps_air*C_air*r)

to become

np.power(10, -eps_air*C_air*r)

Boolean comparison

if no_abs == True should simply be if no_abs.

Enums

model should capture RAD, LSI etc. in an Enum.

set_reactor return value

You're forming a list, which is a little odd. Typically, multi-return values are implicit tuples, i.e.

return stirred_reactor, reactor_network

Spelling

reasults -> results

Redundant return

return None

can be deleted; that's what happens by default.

Dictionary update

plt.rcParams['axes.labelsize'] = 18
plt.rcParams['xtick.labelsize'] = 14
plt.rcParams['ytick.labelsize'] = 14
plt.rcParams['figure.autolayout'] = True

should use

plt.rcParams.update({ ...

Graph bounds

The lower bound for your time axis should be 1e-3 or maybe 1e-2. Remove the long, flat part to focus on the interesting bits.

\$\endgroup\$
2
  • \$\begingroup\$ Thanks got the feedback! I'll apply all the changes but I just wanted some clarification on the Enums. You are suggesting to make a model a class and then use the Enum package to define the different models within the class, right? I don't have much experience with OOP programming, could you give me a pointer to a good resource that I can use as an example? \$\endgroup\$
    – Marc
    Jul 6, 2020 at 8:41
  • \$\begingroup\$ docs.python.org/3/library/enum.html has good examples and is the only "authoritative source" \$\endgroup\$
    – Reinderien
    Jul 6, 2020 at 12:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.