From 2420d97390bc8e6e6071b275c98f3ee13cc27cc2 Mon Sep 17 00:00:00 2001 From: Ga-Yeong Kim <47093338+GaYeongKim@users.noreply.github.com> Date: Mon, 22 Apr 2024 17:44:23 -0500 Subject: [PATCH 01/32] added asm2d components to pm2 components finished working on rho --- qsdsan/data/process_data/_pm2_asm2d.tsv | 39 ++ qsdsan/processes/_pm2_asm2d.py | 793 ++++++++++++++++++++++++ 2 files changed, 832 insertions(+) create mode 100644 qsdsan/data/process_data/_pm2_asm2d.tsv create mode 100644 qsdsan/processes/_pm2_asm2d.py diff --git a/qsdsan/data/process_data/_pm2_asm2d.tsv b/qsdsan/data/process_data/_pm2_asm2d.tsv new file mode 100644 index 00000000..7bc0ff7f --- /dev/null +++ b/qsdsan/data/process_data/_pm2_asm2d.tsv @@ -0,0 +1,39 @@ + X_CHL X_ALG X_CH X_LI S_CO2 S_A S_F S_O2 S_NH S_NO S_P X_N_ALG X_P_ALG S_N2 S_ALK S_I X_I X_S X_H X_AUT +photoadaptation 1 +ammonium_uptake -1 1 +phosphorus_uptake -1 1 +growth_pho 1 ? 1 ? ? +carbohydrate_storage_pho 1 ? 1 +lipid_storage_pho 1 ? 1 +carbohydrate_growth_pho 1 (-Y_CH_PHO/Y_X_ALG_PHO) ? ? ? ? +lipid_growth_pho 1 (-Y_LI_PHO/Y_X_ALG_PHO) ? ? ? ? +carbohydrate_maintenance_pho -1 ? -1 +lipid_maintenance_pho -1 ? -1 +endogenous_respiration_pho -1 ? -1 ? ? +growth_ace 1 ? (-1)/Y_X_ALG_HET_ACE ? ? ? +carbohydrate_storage_ace 1 ? (-1)/Y_CH_ND_HET_ACE ? +lipid_storage_ace 1 ? (-1)/Y_LI_ND_HET_ACE ? +carbohydrate_growth_ace 1 (-Y_CH_NR_HET_ACE/Y_X_ALG_HET_ACE) ? ? ? ? +lipid_growth_ace 1 (-Y_LI_NR_HET_ACE/Y_X_ALG_HET_ACE) ? ? ? ? +carbohydrate_maintenance_ace -1 ? -1 +lipid_maintenance_ace -1 ? -1 +endogenous_respiration_ace -1 ? -1 ? ? +growth_glu 1 ? (-1)/Y_X_ALG_HET_GLU ? ? ? +carbohydrate_storage_glu 1 ? (-1)/Y_CH_ND_HET_GLU ? +lipid_storage_glu 1 ? (-1)/Y_LI_ND_HET_GLU ? +carbohydrate_growth_glu 1 (-Y_CH_NR_HET_GLU/Y_X_ALG_HET_GLU) ? ? ? ? +lipid_growth_glu 1 (-Y_LI_NR_HET_GLU/Y_X_ALG_HET_GLU) ? ? ? ? +carbohydrate_maintenance_glu -1 ? -1 +lipid_maintenance_glu -1 ? -1 +endogenous_respiration_glu -1 ? -1 ? ? +aero_hydrolysis 1-f_SI ? ? ? f_SI -1 +anox_hydrolysis 1-f_SI ? ? ? f_SI -1 +anae_hydrolysis 1-f_SI ? ? ? f_SI -1 +hetero_growth_S_F (-1)/Y_H 1-1/Y_H ? ? ? 1 +hetero_growth_S_A (-1)/Y_H 1-1/Y_H ? ? ? 1 +denitri_S_F (-1)/Y_H ? (Y_H-1)/(20/7*Y_H) ? (1-Y_H)/(20/7*Y_H) ? 1 +denitri_S_A (-1)/Y_H ? (Y_H-1)/(20/7*Y_H) ? (1-Y_H)/(20/7*Y_H) ? 1 +ferment 1 -1 ? ? ? +hetero_lysis ? ? ? f_XI_H 1-f_XI_H -1 +auto_aero_growth (Y_A-32/7)/Y_A ? 1/Y_A ? ? 1 +auto_lysis ? ? ? f_XI_AUT 1-f_XI_AUT -1 diff --git a/qsdsan/processes/_pm2_asm2d.py b/qsdsan/processes/_pm2_asm2d.py new file mode 100644 index 00000000..08a0a294 --- /dev/null +++ b/qsdsan/processes/_pm2_asm2d.py @@ -0,0 +1,793 @@ +# -*- coding: utf-8 -*- +''' +QSDsan: Quantitative Sustainable Design for sanitation and resource recovery systems + +This module is developed by: + Ga-Yeong Kim + Joy Zhang + +This module is under the University of Illinois/NCSA Open Source License. +Please refer to https://github.com/QSD-Group/QSDsan/blob/main/LICENSE.txt +for license details. +''' +# from thermosteam.utils import chemicals_user +from thermosteam import settings +from qsdsan import Component, Components, Process, Processes, CompiledProcesses +from qsdsan.utils import ospath, data_path +import numpy as np + +__all__ = ('create_pm2_asm2d_cmps', 'PM2ASM2d') + +_path = ospath.join(data_path, 'process_data/_pm2_asm2d.tsv') +# _load_components = settings.get_default_chemicals + +#%% +# ============================================================================= +# PM2ASM2d-specific components +# ============================================================================= + +def create_pm2_asm2d_cmps(set_thermo=True): + cmps = Components.load_default() + + # X_CHL (g Chl/m^3) + X_CHL = Component(ID = 'X_CHL', + formula = 'C55H72MgN4O5', + description = 'Chlorophyll content of cells', + particle_size = 'Particulate', + degradability = 'Slowly', + organic = True) + + # X_ALG (g COD/m^3) + X_ALG = cmps.X_OHO.copy('X_ALG') + X_ALG.description = 'Concentration of carbon-accumulating mixotrophic organisms' + X_ALG.formula = 'CH1.8O0.5N0.2P0.018' + X_ALG.f_BOD5_COD = X_ALG.f_uBOD_COD = None + X_ALG.f_Vmass_Totmass = 0.89 + + # X_CH (g COD/m^3) + X_CH = cmps.X_GAO_Gly.copy('X_CH') + X_CH.description = 'Concentration of stored carbohydrates' + X_CH.formula = 'CH2O' + X_CH.f_BOD5_COD = X_CH.f_uBOD_COD = None + + # X_LI (g COD/m^3) + X_LI = cmps.X_GAO_Gly.copy('X_LI') + X_LI.description = 'Concentration of stored lipids' + X_LI.formula = 'CH1.92O0.118' + X_LI.f_BOD5_COD = X_LI.f_uBOD_COD = None + + # S_CO2 (g CO2/m^3) + S_CO2 = Component.from_chemical(ID = 'S_CO2', + chemical = 'CO2', + description = 'Soluble carbon dioxide', + particle_size = 'Soluble', + degradability = 'Undegradable', + organic = False) + + # S_A (g COD/m^3) + S_A = cmps.S_Ac.copy('S_A') + S_A.description = 'Concentration of extracellular dissolved organic carbon (acetate)' + + # S_F (g COD/m^3) + S_F = Component.from_chemical(ID = 'S_F', + chemical = 'glucose', + description = 'Concentration of extracellular dissolved organic carbon (glucose)', + measured_as = 'COD', + particle_size = 'Soluble', + degradability = 'Readily', + organic = True) + + # S_O2 (g O2/m^3) + S_O2 = cmps.S_O2.copy('S_O2') + S_O2.description = ('Concentration of dissolved oxygen') + + # S_NH (g N/m^3) + S_NH = cmps.S_NH4.copy('S_NH') + S_NH.description = ('Concentration of dissolved ammonium') + + # S_NO (g N/m^3) + S_NO = cmps.S_NO3.copy('S_NO') + S_NO.description = ('Concentration of dissolved nitrate/nitrite') + + # S_P (g P/m^3) + S_P = cmps.S_PO4.copy('S_P') + S_P.description = ('Concentration of dissolved phosphorus') + + # X_N_ALG (g N/m^3) + X_N_ALG = cmps.X_B_Subst.copy('X_N_ALG') + X_N_ALG.description = 'Concentration of algal cell-associated nitrogen' + X_N_ALG.measured_as = 'N' + X_N_ALG.i_C = X_N_ALG.i_P = X_N_ALG.i_COD = X_N_ALG.f_BOD5_COD = X_N_ALG.f_uBOD_COD = X_N_ALG.f_Vmass_Totmass = 0 + X_N_ALG.i_mass = 1 + + # X_P_ALG (g P/m^3) + X_P_ALG = cmps.X_B_Subst.copy('X_P_ALG') + X_P_ALG.description = 'Concentration of algal cell-associated phosphorus' + X_P_ALG.measured_as = 'P' + X_P_ALG.i_C = X_P_ALG.i_N = X_P_ALG.i_COD = X_P_ALG.f_BOD5_COD = X_P_ALG.f_uBOD_COD = X_P_ALG.f_Vmass_Totmass = 0 + X_P_ALG.i_mass = 1 + + '''added from asm2d''' + # S_N2 (g N/m^3) + S_N2 = cmps.S_N2.copy('S_N2') + S_N2.description = ('Concentration of dinitrogen') + + # S_ALK (g C/m^3) + S_ALK = cmps.S_CO3.copy('S_ALK') # measured as g C, not as mole HCO3- + S_ALK.description = ('Concentration of alkalinity') + + # S_I (g COD/m^3) + S_I = cmps.S_U_E.copy('S_I') + S_I.description = ('Concentration of inert soluble organic material') + + # X_I (g COD/m^3) + X_I = cmps.X_U_OHO_E.copy('X_I') + X_I.description = ('Concentration of inert particulate organic material') + + # X_S (g COD/m^3) + X_S = cmps.X_B_Subst.copy('X_S') + X_S.description = ('Concentration of slowly biodegradable substrates') + + # X_H (g COD/m^3) + X_H = cmps.X_OHO.copy('X_H') + X_H.description = ('Concentration of heterotrophic organisms (including denitrifer)') + + # X_AUT (g COD/m^3) + X_AUT = cmps.X_AOO.copy('X_AUT') + X_AUT.description = ('Concentration of nitrifying organisms') + + S_I.i_N = 0.01 + # S_F.i_N = 0.03 + X_I.i_N = 0.02 + X_S.i_N = 0.04 + X_H.i_N = X_AUT.i_N = 0.07 + + S_I.i_P = 0.00 + # S_F.i_P = 0.01 + X_I.i_P = 0.01 + X_S.i_P = 0.01 + X_H.i_P = X_AUT.i_P = 0.02 + + X_I.i_mass = 0.75 + X_S.i_mass = 0.75 + X_H.i_mass = X_AUT.i_mass = 0.9 + + cmps_pm2_asm2d = Components([X_CHL, X_ALG, X_CH, X_LI, S_CO2, S_A, S_F, + S_O2, S_NH, S_NO, S_P, X_N_ALG, X_P_ALG, + S_N2, S_ALK, S_I, X_I, X_S, X_H, X_AUT, cmps.H2O]) + + cmps_pm2_asm2d.default_compile() + + if set_thermo: settings.set_thermo(cmps_pm2_asm2d) + return cmps_pm2_asm2d + +# create_pm2_asm2d_cmps() + +#%% +# ============================================================================= +# kinetic rate functions +# ============================================================================= + +# Calculation of ratio +def ratio(numerator, denominator, minimum, maximum): + return min(max(minimum, numerator / denominator), maximum) + +# Calculation of 'I_0' (for initial sensitivity analysis using calculated I_0) +def calc_irrad(t): + ''' + :param t: time [days] + :return: I_0, calculated irradiance [uE/m^2/s] + + -Assumes 14 hours of daylight + ''' + daylight_hours = 14.0 # hours + start_time = (12.0 - daylight_hours / 2) / 24.0 # days + end_time = (12.0 + daylight_hours / 2) / 24.0 # days + if t-np.floor(t) < start_time or t-np.floor(t) > end_time: + return 0 + else: + return 400.0 * (np.sin(2 * np.pi * (((t - np.floor(t)) - 5 / 24) / (14 / 24)) - np.pi / 2) + 1) / 2 + +# Calculation of 'I' from 'I_0' (Beer-Lambert) +def attenuation(light, X_TSS, a_c, b_reactor): + ''' + :param light: I_0, calculated irradiance from 'calc_irrad' method (for sensitivity analysis) or + photosynthetically active radiation (PAR) imported from input excel file (for calibration & validation) [uE/m^2/s] + :param X_TSS: total biomass concentration (X_ALG + X_CH + X_LI) * i_mass [g TSS/m^3] + :param a_c: PAR absorption coefficient on a TSS (total suspended solids) basis [m^2/g TSS] + :parma b_reactor: thickness of reactor along light path [m] + :return: I, depth-averaged irradiance [uE/m^2/s] + ''' + if X_TSS > 0: + i_avg = (light * (1 - np.exp(-a_c * X_TSS * b_reactor))) / (a_c * X_TSS * b_reactor) + return min(i_avg, light) + else: + return light + +# Calculation of 'f_I' from 'I' (Eilers & Peeters) +def irrad_response(i_avg, X_CHL, X_carbon, I_n, I_opt): + ''' + :param i_avg: I, depth-averaged irradiance (calculated from 'attenuation' method) [uE/m^2/s] + :param X_CHL: chlorophyll content of cells [g Chl/m^3] + :param X_carbon: carbon content of cells (X_ALG + X_CH + X_LI) * i_C [g C/m^3] + :param I_n: maximum incident PAR irradiance (“irradiance at noon”) [uE/m^2/s] + :param I_opt: optimal irradiance [uE/m^2/s] + :return: f_I, irradiance response function [unitless] + ''' + if X_carbon > 0: + f_I = i_avg / (i_avg + I_n * (0.25 - (5 * X_CHL/X_carbon)) * ((i_avg ** 2 / I_opt ** 2) - (2 * i_avg / I_opt) + 1)) + return min(1, max(0, f_I)) + else: + return 0 + +# Droop model +def droop(quota, subsistence_quota, exponent): + ''' + :param quota: Q_N or Q_P [g N or g P/g COD] + :param subsistence_quota: Q_N_min or Q_P_min [g N or g P/g COD] + :param exponent: exponent to allow for more rapid transitions from growth to storage (see Guest et al., 2013) [unitless] + :return: rate [unitless] + ''' + return 1 - (subsistence_quota / quota) ** exponent + +# Monod model +def monod(substrate, half_sat_const, exponent): + ''' + :param substrate: S_NH, S_NO or S_P [g N or g P/m^3] + :param half_sat_const: K_N or K_P [g N or g P/m^3] + :param exponent: exponent to allow for more rapid transitions from growth to storage (see Guest et al., 2013) [unitless] + :return: rate [unitless] + ''' + return (substrate / (half_sat_const + substrate)) ** exponent + +# Temperature model (Arrhenius) +def temperature(temp, arr_a, arr_e): + ''' + :param temp: temperature (will be imported from input excel file) [K] + :param arr_a: arrhenius constant (A) (Goldman et al., 1974) [unitless] + :param arr_e: arrhenius exponential constant (E/R) (Goldman et al., 1974) [K] + :return: temperature component of overall growth equation [unitless] + ''' + return arr_a * np.exp(-arr_e / temp) # Used equation from Goldman et al., 1974 + +# Photoadaptation (_p1) +def photoadaptation(i_avg, X_CHL, X_carbon, I_n, k_gamma): + ''' + :param i_avg: I, depth-averaged irradiance (calculated from 'attenuation' method) [uE/m^2/s] + :param X_CHL: chlorophyll content of cells [g Chl/m^3] + :param X_carbon: carbon content of cells (X_ALG + X_CH + X_LI) * i_C [g C/m^3] + :param I_n: maximum incident PAR irradiance (“irradiance at noon”) [uE/m^2/s] + :param k_gamma: photoadaptation coefficient [unitless] + :return: photoadaptation rate [g Chl/m^3/d] + ''' + if X_carbon > 0: + return 24 * ((0.2 * i_avg / I_n) / (k_gamma + (i_avg / I_n))) *\ + (0.01 + 0.03 * ((np.log(i_avg / I_n + 0.005)) / (np.log(0.01))) - X_CHL/X_carbon) * X_carbon + else: return 0 + +# Nutrients uptake (_p2, _p3, _p4, _p5, _p6) +def nutrient_uptake(X_ALG, quota, substrate, uptake_rate, half_sat_const, maximum_quota, subsistence_quota): + ''' + :param X_ALG: algae biomass concentration (i.e., no storage products) [g COD/m^3] + :param quota: Q_N or Q_P [g N or g P/g COD] + :param substrate: S_NH, S_NO or S_P [g N or g P/m^3] + :param uptake_rate: V_NH, V_NO or V_P [g N or g P/g COD/d] + :param half_sat_const: K_N or K_P [g N or g P/m^3] + :param maximum_quota: Q_N_max or Q_P_max [g N or g P/g COD] + :param subsistence_quota: Q_N_min or Q_P_min [g N or g P/g COD] + :return: nutrient uptake rate [g N or g P/m^3/d] + ''' + return uptake_rate * monod(substrate, half_sat_const, 1) * X_ALG * \ + ((maximum_quota - quota) / (maximum_quota - subsistence_quota)) ** 0.01 + +# Maximum total photoautotrophic or heterotrophic-acetate or heterotrophic-glucose growth rate (_p7, _p10, _p11, _p15, _p18, _p19, _p23, _p26, _p27) +def max_total_growth(X_ALG, mu_max, f_np, f_temp): + ''' + :param X_ALG: algae biomass concentration (i.e., no storage products) [g COD/m^3] + :param mu_max: maximum specific growth rate [d^(-1)] + :param f_np: inhibition factor by nitrogen or phosphorus (between 0 and 1) [unitless] + :param f_temp: temperature correction factor (between 0 and 1) [unitless] + :return: maximum total growth rate for a particular mechanism, + without considering carbon source or light inhibition (= product of shared terms in growth-related equations) [g COD/m^3/d] + ''' + return mu_max * f_np * X_ALG * f_temp + +# Split the total growth rate between three processes (_p7, _p10, _p11, _p15, _p18, _p19, _p23, _p26, _p27) +def growth_split(f_I, f_CH, f_LI, rho, Y_CH, Y_LI, K_STO): + ''' + :param f_I: irradiance response function (calculated from 'irrad_response' method) [unitless] + :param f_CH: ratio of stored carbohydrates to cells (X_CH / X_ALG) [g COD/g COD] + :param f_LI: ratio of stored lipids to cells (X_LI / X_ALG) [g COD/g COD] + :param rho: carbohydrate relative preference factor (calibrated in Guest et al., 2013) [unitless] + :param Y_CH: yield of storage carbohydrates (as polyglucose, PG), Y_CH_PHO, Y_CH_NR_HET_ACE, or Y_CH_NR_HET_GLU [g COD/g COD] + :param Y_LI: yield of storage lipids (as triacylglycerol, TAG), Y_LI_PHO, Y_LI_NR_HET_ACE, or Y_LI_NR_HET_GLU [g COD/g COD] + :param K_STO: half-saturation constant for stored organic carbon (calibrated in Guest et al., 2013) [g COD/g COD] + :return: splits the total growth rate between three processes, + growth, growth on stored carbohydrates, and growth on stored lipid (= process-specific terms) [unitless] + ''' + numerators = np.asarray([K_STO * (1 - f_I), rho * f_CH, f_LI * Y_CH / Y_LI]) + return numerators/(sum(numerators)) + +# Part of storage equations (_p8, _p9, _p16, _p17, _p24, _p25) +def storage_saturation(f, f_max, beta): + ''' + :param f: f_CH or f_LI [g COD/g COD] + :param f_max: f_CH_max or f_LI_max [g COD/g COD] + :param beta: beta_1 or beta_2 [unitless] + :return: part of storage equations [unitless] + ''' + return 1 - (f / f_max) ** beta + +# Maximum total photoautotrophic or heterotrophic-acetate or heterotrophic-glucose maintenance rate (_p12, _p13, _p14, _p20, _p21, _p22, _p28, _p29, _p30) +def max_total_maintenance(X_ALG, m_ATP): + ''' + :param X_ALG: algae biomass concentration (i.e., no storage products) [g COD/m^3] + :param m_ATP: specific maintenance rate [g ATP/g COD/d] + :return: maximum total maintenance rate for a particular mechanism + (= product of shared terms in maintenance-related equations) [g COD/m^3/d] + ''' + return m_ATP * X_ALG + +# Split the total maintenance rate between three processes (_p12, _p13, _p14, _p20, _p21, _p22, _p28, _p29, _p30) +def maintenance_split(f_CH, f_LI, rho, Y_CH, Y_LI, Y_X_ALG, Y_ATP, K_STO): + ''' + :param f_CH: ratio of stored carbohydrates to cells (X_CH / X_ALG) [g COD/g COD] + :param f_LI: ratio of stored lipids to cells (X_LI / X_ALG) [g COD/g COD] + :param rho: carbohydrate relative preference factor (calibrated in Guest et al., 2013) [unitless] + :param Y_CH: yield of storage carbohydrates (as polyglucose, PG), Y_CH_PHO, Y_CH_NR_HET_ACE, or Y_CH_NR_HET_GLU [g COD/g COD] + :param Y_LI: yield of storage lipids (as triacylglycerol, TAG), Y_LI_PHO, Y_LI_NR_HET_ACE, or Y_LI_NR_HET_GLU [g COD/g COD] + :param Y_X_ALG: yield of carbon-accumulating phototrophic organisms, Y_X_ALG_PHO, Y_X_ALG_HET_ACE, or Y_X_ALG_HET_GLU [g COD/g COD] + :param Y_ATP: yield of ATP, Y_ATP_PHO, Y_ATP_HET_ACE, or Y_ATP_HET_GLU [g ATP/g COD] + :param K_STO: half-saturation constant for stored organic carbon (calibrated in Guest et al., 2013) [g COD/g COD] + :return: splits the total maintenance rate between three processes, + stored carbohydrate degradation, stored lipid degradation, and endogenous respiration (= process-specific terms) [unitless] + ''' + numerators = np.asarray([rho * f_CH, f_LI * Y_CH / Y_LI, K_STO]) + yield_ratios = np.asarray([Y_CH, Y_LI, Y_X_ALG]) / Y_ATP + return numerators/(sum(numerators)) * yield_ratios + +# Storage of carbohydrate/lipid (_p8, _p9, _p16, _p17, _p24, _p25) +def storage(X_ALG, f_np, response, saturation, storage_rate): + ''' + :param X_ALG: algae biomass concentration (i.e., no storage products) [g COD/m^3] + :param f_np: inhibition factor by nitrogen or phosphorus (between 0 and 1) [unitless] + :param response: f_I (irradiance response function, calculated from 'irrad_response' method), acetate_response (monod(S_A, K_A, 1)), or glucose_response (monod(S_F, K_F, 1)) [unitless] + :param saturation: 1 - (f / f_max) ** beta (calculated from 'storage_saturation' method) [unitless] + :param storage_rate: q_CH or q_LI [g COD/g COD/d] + :return: storage rate [g COD/m^3/d] + ''' + return storage_rate * saturation * (1 - f_np) * response * X_ALG + +def rhos_pm2_asm2d(state_arr, params): + + # extract values of state variables + c_arr = state_arr[:21] + temp = state_arr[22] + light = state_arr[23] # imported from input file assumed + + # Q = state_arr[14] # Flow rate + # t = state_arr[15] # time + + # light = calc_irrad(t) # when to use calculated light (I_0) + + # X_CHL, X_ALG, X_CH, X_LI, S_CO2, S_A, S_F, S_O2, S_NH, S_NO, S_P, X_N_ALG, X_P_ALG, H2O = c_arr + X_CHL, X_ALG, X_CH, X_LI, S_CO2, S_A, S_F, S_O2, S_NH, S_NO, S_P, X_N_ALG, X_P_ALG, S_N2, S_ALK, S_I, X_I, X_S, X_H, X_AUT, H2O = c_arr + + # extract values of parameters + cmps = params['cmps'] + a_c = params['a_c'] + I_n = params['I_n'] + arr_a = params['arr_a'] + arr_e = params['arr_e'] + beta_1 = params['beta_1'] + beta_2 = params['beta_2'] + b_reactor = params['b_reactor'] + I_opt = params['I_opt'] + k_gamma = params['k_gamma'] + K_N = params['K_N'] + K_P = params['K_P'] + K_A = params['K_A'] + K_F = params['K_F'] + rho = params['rho'] + K_STO = params['K_STO'] + f_CH_max = params['f_CH_max'] + f_LI_max = params['f_LI_max'] + m_ATP = params['m_ATP'] + mu_max = params['mu_max'] + q_CH = params['q_CH'] + q_LI = params['q_LI'] + Q_N_max = params['Q_N_max'] + Q_N_min = params['Q_N_min'] + Q_P_max = params['Q_P_max'] + Q_P_min = params['Q_P_min'] + V_NH = params['V_NH'] + V_NO = params['V_NO'] + V_P = params['V_P'] + exponent = params['exponent'] + Y_ATP_PHO = params['Y_ATP_PHO'] + Y_CH_PHO = params['Y_CH_PHO'] + Y_LI_PHO = params['Y_LI_PHO'] + Y_X_ALG_PHO = params['Y_X_ALG_PHO'] + Y_ATP_HET_ACE = params['Y_ATP_HET_ACE'] + Y_CH_NR_HET_ACE = params['Y_CH_NR_HET_ACE'] + Y_LI_NR_HET_ACE = params['Y_LI_NR_HET_ACE'] + Y_X_ALG_HET_ACE = params['Y_X_ALG_HET_ACE'] + Y_ATP_HET_GLU = params['Y_ATP_HET_GLU'] + Y_CH_NR_HET_GLU = params['Y_CH_NR_HET_GLU'] + Y_LI_NR_HET_GLU = params['Y_LI_NR_HET_GLU'] + Y_X_ALG_HET_GLU = params['Y_X_ALG_HET_GLU'] + n_dark = params['n_dark'] + + '''added from asm2d''' + f_SI = params['f_SI'] + Y_H = params['Y_H'] + f_XI_H = params['f_XI_H'] + Y_A = params['Y_A'] + f_XI_AUT = params['f_XI_AUT'] + K_h = params['K_h'] + eta_NO3 = params['eta_NO3'] + eta_fe = params['eta_fe'] + K_O2 = params['K_O2'] + K_NO3 = params['K_NO3'] + K_X = params['K_X'] + mu_H = params['mu_H'] + q_fe = params['q_fe'] + eta_NO3_H = params['eta_NO3_H'] + b_H = params['b_H'] + K_O2_H = params['K_O2_H'] + K_F_H = params['K_F_H'] # K_F overlaps with PM2 -> change into K_F_H + K_fe = params['K_fe'] + K_A_H = params['K_A_H'] + K_NO3_H = params['K_NO3_H'] + K_NH4_H = params['K_NH4_H'] + K_P_H = params['K_P_H'] + K_ALK_H = params['K_ALK_H'] + mu_AUT = params['mu_AUT'] + b_AUT = params['b_AUT'] + K_O2_AUT = params['K_O2_AUT'] + K_NH4_AUT = params['K_NH4_AUT'] + K_ALK_AUT = params['K_ALK_AUT'] + K_P_AUT = params['K_P_AUT'] + +# intermediate variables + f_CH = ratio(X_CH, X_ALG, 0, f_CH_max) + f_LI = ratio(X_LI, X_ALG, 0, f_LI_max) + + # Q_N = ratio(X_N_ALG, X_ALG, Q_N_min, Q_N_max) + # Q_P = ratio(X_P_ALG, X_ALG, Q_P_min, Q_P_max) + + alg_iN, alg_iP = cmps.X_ALG.i_N, cmps.X_ALG.i_P + Q_N = ratio(X_N_ALG+X_ALG*alg_iN, X_ALG, Q_N_min, Q_N_max) + Q_P = ratio(X_P_ALG+X_ALG*alg_iP, X_ALG, Q_P_min, Q_P_max) + + idx = cmps.indices(['X_ALG', 'X_CH', 'X_LI']) + X_bio = np.array([X_ALG, X_CH, X_LI]) + X_TSS = sum(X_bio * cmps.i_mass[idx]) + X_carbon = sum(X_bio * cmps.i_C[idx]) + + i_avg = attenuation(light, X_TSS, a_c, b_reactor) + f_I = irrad_response(i_avg, X_CHL, X_carbon, I_n, I_opt) + dark_response = max(f_I, n_dark) + acetate_response = monod(S_A, K_A, 1) + glucose_response = monod(S_F, K_F, 1) + + f_np = min(droop(Q_N, Q_N_min, exponent), droop(Q_P, Q_P_min, exponent)) + f_temp = temperature(temp, arr_a, arr_e) + + f_sat_CH = storage_saturation(f_CH, f_CH_max, beta_1) + f_sat_LI = storage_saturation(f_LI, f_LI_max, beta_2) + + max_total_growth_rho = max_total_growth(X_ALG, mu_max, f_np, f_temp) + max_maintenance_rho = max_total_maintenance(X_ALG, m_ATP) + # light = calc_irrad(t) + + # calculate kinetic rate values + rhos = np.empty(30) + + rhos[0] = photoadaptation(i_avg, X_CHL, X_carbon, I_n, k_gamma) + + rhos[1] = nutrient_uptake(X_ALG, Q_N, S_NH, V_NH, K_N, Q_N_max, Q_N_min) + rhos[[2,3,4]] = nutrient_uptake(X_ALG, Q_N, S_NO, V_NO, K_N, Q_N_max, Q_N_min) * (K_N/(K_N + S_NH)) + + rhos[5] = nutrient_uptake(X_ALG, Q_P, S_P, V_P, K_P, Q_P_max, Q_P_min) + + rhos[[6,9,10]] = max_total_growth_rho \ + * growth_split(f_I, f_CH, f_LI, rho, Y_CH_PHO, Y_LI_PHO, K_STO) + rhos[6] *= f_I + rhos[[9,10]] *= dark_response + + rhos[[14,17,18]] = max_total_growth_rho \ + * acetate_response \ + * growth_split(f_I, f_CH, f_LI, rho, Y_CH_NR_HET_ACE, Y_LI_NR_HET_ACE, K_STO) + + rhos[[22,25,26]] = max_total_growth_rho \ + * glucose_response \ + * growth_split(f_I, f_CH, f_LI, rho, Y_CH_NR_HET_GLU, Y_LI_NR_HET_GLU, K_STO) + + rhos[[11,12,13]] = max_maintenance_rho \ + * maintenance_split(f_CH, f_LI, rho, Y_CH_PHO, Y_LI_PHO, Y_X_ALG_PHO, Y_ATP_PHO, K_STO) + + rhos[[19,20,21]] = max_maintenance_rho \ + * maintenance_split(f_CH, f_LI, rho, Y_CH_NR_HET_ACE, Y_LI_NR_HET_ACE, Y_X_ALG_HET_ACE, Y_ATP_HET_ACE, K_STO) + + rhos[[27,28,29]] = max_maintenance_rho \ + * maintenance_split(f_CH, f_LI, rho, Y_CH_NR_HET_GLU, Y_LI_NR_HET_GLU, Y_X_ALG_HET_GLU, Y_ATP_HET_GLU, K_STO) + + rhos[7] = storage(X_ALG, f_np, f_I, f_sat_CH, q_CH) + rhos[8] = storage(X_ALG, f_np, f_I, f_sat_LI, q_LI) * (f_CH / f_CH_max) + rhos[15] = storage(X_ALG, f_np, acetate_response, f_sat_CH, q_CH) + rhos[16] = storage(X_ALG, f_np, acetate_response, f_sat_LI, q_LI) * (f_CH / f_CH_max) + rhos[23] = storage(X_ALG, f_np, glucose_response, f_sat_CH, q_CH) + rhos[24] = storage(X_ALG, f_np, glucose_response, f_sat_LI, q_LI) * (f_CH / f_CH_max) + + return rhos + +#%% +# ============================================================================= +# PM2 class +# ============================================================================= + +class PM2ASM2d(CompiledProcesses): + ''' + Parameters + ---------- + components: class:`CompiledComponents`, optional + Components corresponding to each entry in the stoichiometry array, + defaults to thermosteam.settings.chemicals. + a_c : float, optional + PAR absorption coefficient on a TSS (total suspended solids) basis, in [m^2/g TSS]. + The default is 0.049. + I_n : float, optional + Maximum incident PAR irradiance (“irradiance at noon”), in [uE/m^2/s]. + The default is 250. + arr_a : float, optional + Arrhenius constant (A), in [unitless]. + The default is 1.8 * 10**10. + arr_e : float, optional + Arrhenius exponential constant (E/R), in [K]. + The default is 6842. + beta_1 : float, optional + Power coefficient for carbohydrate storage inhibition, in [unitless]. + The default is 2.90. + beta_2 : float, optional + Power coefficient for lipid storage inhibition, in [unitless]. + The default is 3.50. + b_reactor : float, optional + Thickness of reactor along light path, in [m]. + The default is 0.03. + I_opt : float, optional + Optimal irradiance, in [uE/m^2/s]. + The default is 300. + k_gamma : float, optional + Photoadaptation coefficient, in [unitless]. + The default is 0.00001. + K_N : float, optional + Nitrogen half-saturation constant, in [g N/m^3]. + The default is 0.1. + K_P : float, optional + Phosphorus half-saturation constant, in [g P/m^3]. + The default is 1.0. + K_A : float, optional + Organic carbon half-saturation constant (acetate) (Wagner, 2016), in [g COD/m^3]. + The default is 6.3. + K_F : float, optional + Organic carbon half-saturation constant (glucose); assumes K_A = K_F, in [g COD/m^3]. + The default is 6.3. + rho : float, optional + Carbohydrate relative preference factor (calibrated in Guest et al., 2013), in [unitless]. + The default is 1.186. + K_STO : float, optional + Half-saturation constant for stored organic carbon (calibrated in Guest et al., 2013), in [g COD/g COD]. + The default is 1.566. + f_CH_max : float, optional + Maximum achievable ratio of stored carbohydrates to functional cells, in [g COD/g COD]. + The default is 0.819. + f_LI_max : float, optional + Maximum achievable ratio of stored lipids to functional cells, in [g COD/g COD]. + The default is 3.249. + m_ATP : float, optional + Specific maintenance rate, in [g ATP/g COD/d]. + The default is 15.835. + mu_max : float, optional + Maximum specific growth rate, in [d^(-1)]. + The default is 1.969. + q_CH : float, optional + Maximum specific carbohydrate storage rate, in [g COD/g COD/d]. + The default is 0.594. + q_LI : float, optional + Maximum specific lipid storage rate, in [g COD/g COD/d]. + The default is 0.910. + Q_N_max : float, optional + Maximum nitrogen quota, in [g N/g COD]. + The default is 0.417. + Q_N_min : float, optional + Nitrogen subsistence quota, in [g N/g COD]. + The default is 0.082. + Q_P_max : float, optional + Maximum phosphorus quota, in [g P/g COD]. + The default is 0.092. + Q_P_min : float, optional + Phosphorus subsistence quota; assumes N:P ratio of 5:1, in [g P/g COD]. + The default is 0.0163. + V_NH : float, optional + Maximum specific ammonium uptake rate (calibrated in Guest et al., 2013), in [g N/g COD/d]. + The default is 0.254. + V_NO : float, optional + Maximum specific nitrate uptake rate (calibrated in Guest et al., 2013), in [g N/g COD/d]. + The default is 0.254. + V_P : float, optional + Maximum specific phosphorus uptake rate (calibrated in Guest et al., 2013), in [g P/g COD/d]. + The default is 0.016. + exponent : float, optional + Exponent to allow for more rapid transitions from growth to storage (see Guest et al., 2013), in [unitless] + The default is 4. + Y_ATP_PHO : float, optional + Yield of ATP on CO2 fixed to G3P, in [g ATP/g CO2]. + The default is 55.073. + Y_CH_PHO : float, optional + Yield of storage carbohydrate (as polyglucose, PG) on CO2 fixed to G3P, in [g COD/g CO2]. + The default is 0.754. + Y_LI_PHO : float, optional + Yield of storage lipids (as triacylglycerol, TAG) on CO2 fixed to G3P, in [g COD/g CO2]. + The default is 0.901. + Y_X_ALG_PHO : float, optional + Yield of carbon-accumulating phototrophic organisms on CO2 fixed to G3P, in [g COD/g CO2]. + The default is 0.450. + Y_ATP_HET_ACE : float, optional + Yield of ATP on acetate fixed to acetyl-CoA, in [g ATP/g COD]. + The default is 39.623. + Y_CH_NR_HET_ACE : float, optional + Yield of storage carbohydrates (as polyglucose, PG) on acetate fixed to acetyl-CoA under nutrient-replete condition, in [g COD/g COD]. + The default is 0.625. + Y_CH_ND_HET_ACE : float, optional + Yield of storage carbohydrates (as polyglucose, PG) on acetate fixed to acetyl-CoA under nutrient-deplete condition, in [g COD/g COD]. + The default is 0.600. + Y_LI_NR_HET_ACE : float, optional + Yield of storage lipids (as triacylglycerol, TAG) on acetate fixed to acetyl-CoA under nutrient-replete condition, in [g COD/g COD]. + The default is 1.105. + Y_LI_ND_HET_ACE : float, optional + Yield of storage lipids (as triacylglycerol, TAG) on acetate fixed to acetyl-CoA under nutrient-deplete condition, in [g COD/g COD]. + The default is 0.713. + Y_X_ALG_HET_ACE : float, optional + Yield of carbon-accumulating phototrophic organisms on acetate fixed to acetyl-CoA, in [g COD/g COD]. + The default is 0.216. + Y_ATP_HET_GLU : float, optional + Yield of ATP on glucose fixed to G6P, in [g ATP/g COD]. + The default is 58.114. + Y_CH_NR_HET_GLU : float, optional + Yield of storage carbohydrates (as polyglucose, PG) on glucose fixed to G6P under nutrient-replete condition, in [g COD/g COD]. + The default is 0.917. + Y_CH_ND_HET_GLU : float, optional + Yield of storage carbohydrates (as polyglucose, PG) on glucose fixed to G6P under nutrient-deplete condition, in [g COD/g COD]. + The default is 0.880. + Y_LI_NR_HET_GLU : float, optional + Yield of storage lipids (as triacylglycerol, TAG) on glucose fixed to G6P under nutrient-replete condition, in [g COD/g COD]. + The default is 1.620. + Y_LI_ND_HET_GLU : float, optional + Yield of storage lipids (as triacylglycerol, TAG) on glucose fixed to G6P under nutrient-deplete condition, in [g COD/g COD]. + The default is 1.046. + Y_X_ALG_HET_GLU : float, optional + Yield of carbon-accumulating phototrophic organisms on glucose fixed to G6P, in [g COD/g COD]. + The default is 0.317. + n_dark: float, optional + Dark growth reduction factor, in [unitless] + The default is 0.7. + path : str, optional + Alternative file path for the Petersen matrix. + The default is None. + + Examples + -------- + >>> from qsdsan import processes as pc + >>> cmps = pc.create_pm2_cmps() + >>> pm2 = pc.PM2() + >>> pm2.show() + PM2([photoadaptation, ammonium_uptake, nitrate_uptake_pho, nitrate_uptake_ace, nitrate_uptake_glu, phosphorus_uptake, + growth_pho, carbohydrate_storage_pho, lipid_storage_pho, carbohydrate_growth_pho, lipid_growth_pho, + carbohydrate_maintenance_pho, lipid_maintenance_pho, endogenous_respiration_pho, + growth_ace, carbohydrate_storage_ace, lipid_storage_ace, carbohydrate_growth_ace, lipid_growth_ace, + carbohydrate_maintenance_ace, lipid_maintenance_ace, endogenous_respiration_ace, + growth_glu, carbohydrate_storage_glu, lipid_storage_glu, carbohydrate_growth_glu, lipid_growth_glu, + carbohydrate_maintenance_glu, lipid_maintenance_glu, endogenous_respiration_glu]) + ''' + + _shared_params = ('Y_CH_PHO', 'Y_LI_PHO', 'Y_X_ALG_PHO', + 'Y_CH_NR_HET_ACE', 'Y_LI_NR_HET_ACE', 'Y_X_ALG_HET_ACE', + 'Y_CH_NR_HET_GLU', 'Y_LI_NR_HET_GLU', 'Y_X_ALG_HET_GLU') + + _stoichio_params = ('Y_CH_ND_HET_ACE', 'Y_LI_ND_HET_ACE', 'Y_CH_ND_HET_GLU', 'Y_LI_ND_HET_GLU', + *_shared_params) + + _kinetic_params = ('a_c', 'I_n', 'arr_a', 'arr_e', 'beta_1', 'beta_2', 'b_reactor', 'I_opt', 'k_gamma', + 'K_N', 'K_P', 'K_A', 'K_F', 'rho', 'K_STO', 'f_CH_max', 'f_LI_max', 'm_ATP', 'mu_max', + 'q_CH', 'q_LI', 'Q_N_max', 'Q_N_min', 'Q_P_max', 'Q_P_min', 'V_NH', 'V_NO', 'V_P', 'exponent', + 'Y_ATP_PHO', 'Y_ATP_HET_ACE', 'Y_ATP_HET_GLU', *_shared_params, 'n_dark', 'cmps') + + def __new__(cls, components=None, + a_c=0.049, I_n=250, arr_a=1.8e10, arr_e=6842, beta_1=2.90, beta_2=3.50, b_reactor=0.03, I_opt=300, k_gamma=1e-5, + K_N=0.1, K_P=1.0, K_A=6.3, K_F=6.3, rho=1.186, K_STO=1.566, + f_CH_max=0.819, f_LI_max=3.249, m_ATP=15.835, mu_max=1.969, q_CH=0.594, q_LI=0.910, + Q_N_max=0.417, Q_N_min=0.082, Q_P_max=0.092, Q_P_min=0.0163, V_NH=0.254, V_NO=0.254, V_P=0.016, exponent=4, + Y_ATP_PHO=55.073, Y_CH_PHO=0.754, Y_LI_PHO=0.901, Y_X_ALG_PHO=0.450, + Y_ATP_HET_ACE=39.623, Y_CH_NR_HET_ACE=0.625, Y_CH_ND_HET_ACE=0.600, + Y_LI_NR_HET_ACE=1.105, Y_LI_ND_HET_ACE=0.713, Y_X_ALG_HET_ACE=0.216, + Y_ATP_HET_GLU=58.114, Y_CH_NR_HET_GLU=0.917, Y_CH_ND_HET_GLU=0.880, + Y_LI_NR_HET_GLU=1.620, Y_LI_ND_HET_GLU=1.046, Y_X_ALG_HET_GLU=0.317, n_dark=0.7, + path=None, **kwargs): + + if not path: path = _path + self = Processes.load_from_file(path, + components=components, + conserved_for=('COD', 'C', 'N', 'P'), + parameters=cls._stoichio_params, + compile=False) + + if path == _path: + _p3 = Process('nitrate_uptake_pho', + 'S_NO -> [?]S_O2 + X_N_ALG', + components=components, + ref_component='X_N_ALG', + conserved_for=('COD', 'C')) + + _p4 = Process('nitrate_uptake_ace', + 'S_NO + [?]S_A -> [?]S_CO2 + X_N_ALG', + components=components, + ref_component='X_N_ALG', + conserved_for=('COD', 'C')) + + _p5 = Process('nitrate_uptake_glu', + 'S_NO + [?]S_F -> [?]S_CO2 + X_N_ALG', + components=components, + ref_component='X_N_ALG', + conserved_for=('COD', 'C')) + + self.insert(2, _p3) + self.insert(3, _p4) + self.insert(4, _p5) + + self.compile(to_class=cls) + + self.set_rate_function(rhos_pm2_asm2d) + shared_values = (Y_CH_PHO, Y_LI_PHO, Y_X_ALG_PHO, + Y_CH_NR_HET_ACE, Y_LI_NR_HET_ACE, Y_X_ALG_HET_ACE, + Y_CH_NR_HET_GLU, Y_LI_NR_HET_GLU, Y_X_ALG_HET_GLU) + stoichio_values = (Y_CH_ND_HET_ACE, Y_LI_ND_HET_ACE, Y_CH_ND_HET_GLU, Y_LI_ND_HET_GLU, + *shared_values) + Q_N_min = max(self.Th_Q_N_min, Q_N_min) + Q_P_min = max(self.Th_Q_P_min, Q_P_min) + kinetic_values = (a_c, I_n, arr_a, arr_e, beta_1, beta_2, b_reactor, I_opt, k_gamma, + K_N, K_P, K_A, K_F, rho, K_STO, f_CH_max, f_LI_max, m_ATP, mu_max, + q_CH, q_LI, Q_N_max, Q_N_min, Q_P_max, Q_P_min, V_NH, V_NO, V_P, exponent, + Y_ATP_PHO, Y_ATP_HET_ACE, Y_ATP_HET_GLU, + *shared_values, n_dark, self._components) + + dct = self.__dict__ + dct.update(kwargs) + dct['_parameters'] = dict(zip(cls._stoichio_params, stoichio_values)) + self.rate_function._params = dict(zip(cls._kinetic_params, kinetic_values)) + + return self + + def set_parameters(self, **parameters): + '''Set values to stoichiometric and/or kinetic parameters.''' + stoichio_only = {k:v for k,v in parameters.items() if k in self._stoichio_params} + self._parameters.update(stoichio_only) + if self._stoichio_lambdified is not None: + self.__dict__['_stoichio_lambdified'] = None + if 'Q_N_min' in parameters.keys(): + if parameters['Q_N_min'] < self.Th_Q_N_min: + raise ValueError(f'Value for Q_N_min must not be less than the ' + f'theoretical minimum {self.Th_Q_N_min}') + if 'Q_P_min' in parameters.keys(): + if parameters['Q_P_min'] < self.Th_Q_P_min: + raise ValueError(f'Value for Q_P_min must not be less than the ' + f'theoretical minimum {self.Th_Q_P_min}') + self.rate_function.set_param(**parameters) + + @property + def Th_Q_N_min(self): + return abs(self.stoichiometry.loc['growth_pho', 'X_N_ALG'])*1.001 + + @property + def Th_Q_P_min(self): + return abs(self.stoichiometry.loc['growth_pho', 'X_P_ALG'])*1.001 \ No newline at end of file From d9147984eaec8cd1df1938182398db2f406489ee Mon Sep 17 00:00:00 2001 From: Ga-Yeong Kim Date: Mon, 22 Apr 2024 22:26:31 -0500 Subject: [PATCH 02/32] added rate equations rate equations updated final class underway --- .../{_pm2_asm2d.tsv => _pm2asm2d.tsv} | 0 .../processes/{_pm2_asm2d.py => _pm2asm2d.py} | 404 +++++++++++------- 2 files changed, 243 insertions(+), 161 deletions(-) rename qsdsan/data/process_data/{_pm2_asm2d.tsv => _pm2asm2d.tsv} (100%) rename qsdsan/processes/{_pm2_asm2d.py => _pm2asm2d.py} (79%) diff --git a/qsdsan/data/process_data/_pm2_asm2d.tsv b/qsdsan/data/process_data/_pm2asm2d.tsv similarity index 100% rename from qsdsan/data/process_data/_pm2_asm2d.tsv rename to qsdsan/data/process_data/_pm2asm2d.tsv diff --git a/qsdsan/processes/_pm2_asm2d.py b/qsdsan/processes/_pm2asm2d.py similarity index 79% rename from qsdsan/processes/_pm2_asm2d.py rename to qsdsan/processes/_pm2asm2d.py index 08a0a294..7a0d98b2 100644 --- a/qsdsan/processes/_pm2_asm2d.py +++ b/qsdsan/processes/_pm2asm2d.py @@ -16,97 +16,97 @@ from qsdsan.utils import ospath, data_path import numpy as np -__all__ = ('create_pm2_asm2d_cmps', 'PM2ASM2d') +__all__ = ('create_pm2asm2d_cmps', 'PM2ASM2d') -_path = ospath.join(data_path, 'process_data/_pm2_asm2d.tsv') +_path = ospath.join(data_path, 'process_data/_pm2asm2d.tsv') # _load_components = settings.get_default_chemicals -#%% +#%% # ============================================================================= # PM2ASM2d-specific components # ============================================================================= -def create_pm2_asm2d_cmps(set_thermo=True): +def create_pm2asm2d_cmps(set_thermo=True): cmps = Components.load_default() # X_CHL (g Chl/m^3) X_CHL = Component(ID = 'X_CHL', - formula = 'C55H72MgN4O5', + formula = 'C55H72MgN4O5', description = 'Chlorophyll content of cells', - particle_size = 'Particulate', - degradability = 'Slowly', - organic = True) - + particle_size = 'Particulate', + degradability = 'Slowly', + organic = True) + # X_ALG (g COD/m^3) - X_ALG = cmps.X_OHO.copy('X_ALG') + X_ALG = cmps.X_OHO.copy('X_ALG') X_ALG.description = 'Concentration of carbon-accumulating mixotrophic organisms' - X_ALG.formula = 'CH1.8O0.5N0.2P0.018' + X_ALG.formula = 'CH1.8O0.5N0.2P0.018' X_ALG.f_BOD5_COD = X_ALG.f_uBOD_COD = None - X_ALG.f_Vmass_Totmass = 0.89 - + X_ALG.f_Vmass_Totmass = 0.89 + # X_CH (g COD/m^3) - X_CH = cmps.X_GAO_Gly.copy('X_CH') + X_CH = cmps.X_GAO_Gly.copy('X_CH') X_CH.description = 'Concentration of stored carbohydrates' X_CH.formula = 'CH2O' X_CH.f_BOD5_COD = X_CH.f_uBOD_COD = None - + # X_LI (g COD/m^3) - X_LI = cmps.X_GAO_Gly.copy('X_LI') + X_LI = cmps.X_GAO_Gly.copy('X_LI') X_LI.description = 'Concentration of stored lipids' X_LI.formula = 'CH1.92O0.118' X_LI.f_BOD5_COD = X_LI.f_uBOD_COD = None - + # S_CO2 (g CO2/m^3) - S_CO2 = Component.from_chemical(ID = 'S_CO2', + S_CO2 = Component.from_chemical(ID = 'S_CO2', chemical = 'CO2', description = 'Soluble carbon dioxide', particle_size = 'Soluble', - degradability = 'Undegradable', - organic = False) - + degradability = 'Undegradable', + organic = False) + # S_A (g COD/m^3) - S_A = cmps.S_Ac.copy('S_A') + S_A = cmps.S_Ac.copy('S_A') S_A.description = 'Concentration of extracellular dissolved organic carbon (acetate)' - + # S_F (g COD/m^3) S_F = Component.from_chemical(ID = 'S_F', - chemical = 'glucose', + chemical = 'glucose', description = 'Concentration of extracellular dissolved organic carbon (glucose)', - measured_as = 'COD', + measured_as = 'COD', particle_size = 'Soluble', - degradability = 'Readily', - organic = True) - + degradability = 'Readily', + organic = True) + # S_O2 (g O2/m^3) - S_O2 = cmps.S_O2.copy('S_O2') - S_O2.description = ('Concentration of dissolved oxygen') - + S_O2 = cmps.S_O2.copy('S_O2') + S_O2.description = ('Concentration of dissolved oxygen') + # S_NH (g N/m^3) - S_NH = cmps.S_NH4.copy('S_NH') - S_NH.description = ('Concentration of dissolved ammonium') - + S_NH = cmps.S_NH4.copy('S_NH') + S_NH.description = ('Concentration of dissolved ammonium') + # S_NO (g N/m^3) - S_NO = cmps.S_NO3.copy('S_NO') + S_NO = cmps.S_NO3.copy('S_NO') S_NO.description = ('Concentration of dissolved nitrate/nitrite') - + # S_P (g P/m^3) - S_P = cmps.S_PO4.copy('S_P') + S_P = cmps.S_PO4.copy('S_P') S_P.description = ('Concentration of dissolved phosphorus') - + # X_N_ALG (g N/m^3) X_N_ALG = cmps.X_B_Subst.copy('X_N_ALG') X_N_ALG.description = 'Concentration of algal cell-associated nitrogen' X_N_ALG.measured_as = 'N' X_N_ALG.i_C = X_N_ALG.i_P = X_N_ALG.i_COD = X_N_ALG.f_BOD5_COD = X_N_ALG.f_uBOD_COD = X_N_ALG.f_Vmass_Totmass = 0 X_N_ALG.i_mass = 1 - + # X_P_ALG (g P/m^3) X_P_ALG = cmps.X_B_Subst.copy('X_P_ALG') X_P_ALG.description = 'Concentration of algal cell-associated phosphorus' X_P_ALG.measured_as = 'P' X_P_ALG.i_C = X_P_ALG.i_N = X_P_ALG.i_COD = X_P_ALG.f_BOD5_COD = X_P_ALG.f_uBOD_COD = X_P_ALG.f_Vmass_Totmass = 0 X_P_ALG.i_mass = 1 - + '''added from asm2d''' # S_N2 (g N/m^3) S_N2 = cmps.S_N2.copy('S_N2') @@ -151,18 +151,18 @@ def create_pm2_asm2d_cmps(set_thermo=True): X_I.i_mass = 0.75 X_S.i_mass = 0.75 X_H.i_mass = X_AUT.i_mass = 0.9 - - cmps_pm2_asm2d = Components([X_CHL, X_ALG, X_CH, X_LI, S_CO2, S_A, S_F, - S_O2, S_NH, S_NO, S_P, X_N_ALG, X_P_ALG, + + cmps_pm2asm2d = Components([X_CHL, X_ALG, X_CH, X_LI, S_CO2, S_A, S_F, + S_O2, S_NH, S_NO, S_P, X_N_ALG, X_P_ALG, S_N2, S_ALK, S_I, X_I, X_S, X_H, X_AUT, cmps.H2O]) - - cmps_pm2_asm2d.default_compile() - if set_thermo: settings.set_thermo(cmps_pm2_asm2d) - return cmps_pm2_asm2d + cmps_pm2asm2d.default_compile() + + if set_thermo: settings.set_thermo(cmps_pm2asm2d) + return cmps_pm2asm2d + +# create_pm2asm2d_cmps() -# create_pm2_asm2d_cmps() - #%% # ============================================================================= # kinetic rate functions @@ -191,11 +191,11 @@ def calc_irrad(t): # Calculation of 'I' from 'I_0' (Beer-Lambert) def attenuation(light, X_TSS, a_c, b_reactor): ''' - :param light: I_0, calculated irradiance from 'calc_irrad' method (for sensitivity analysis) or - photosynthetically active radiation (PAR) imported from input excel file (for calibration & validation) [uE/m^2/s] + :param light: I_0, calculated irradiance from 'calc_irrad' method (for sensitivity analysis) or + photosynthetically active radiation (PAR) imported from input excel file (for calibration & validation) [uE/m^2/s] :param X_TSS: total biomass concentration (X_ALG + X_CH + X_LI) * i_mass [g TSS/m^3] - :param a_c: PAR absorption coefficient on a TSS (total suspended solids) basis [m^2/g TSS] - :parma b_reactor: thickness of reactor along light path [m] + :param a_c: PAR absorption coefficient on a TSS (total suspended solids) basis [m^2/g TSS] + :parma b_reactor: thickness of reactor along light path [m] :return: I, depth-averaged irradiance [uE/m^2/s] ''' if X_TSS > 0: @@ -208,9 +208,9 @@ def attenuation(light, X_TSS, a_c, b_reactor): def irrad_response(i_avg, X_CHL, X_carbon, I_n, I_opt): ''' :param i_avg: I, depth-averaged irradiance (calculated from 'attenuation' method) [uE/m^2/s] - :param X_CHL: chlorophyll content of cells [g Chl/m^3] - :param X_carbon: carbon content of cells (X_ALG + X_CH + X_LI) * i_C [g C/m^3] - :param I_n: maximum incident PAR irradiance (“irradiance at noon”) [uE/m^2/s] + :param X_CHL: chlorophyll content of cells [g Chl/m^3] + :param X_carbon: carbon content of cells (X_ALG + X_CH + X_LI) * i_C [g C/m^3] + :param I_n: maximum incident PAR irradiance (“irradiance at noon”) [uE/m^2/s] :param I_opt: optimal irradiance [uE/m^2/s] :return: f_I, irradiance response function [unitless] ''' @@ -227,7 +227,7 @@ def droop(quota, subsistence_quota, exponent): :param subsistence_quota: Q_N_min or Q_P_min [g N or g P/g COD] :param exponent: exponent to allow for more rapid transitions from growth to storage (see Guest et al., 2013) [unitless] :return: rate [unitless] - ''' + ''' return 1 - (subsistence_quota / quota) ** exponent # Monod model @@ -240,25 +240,25 @@ def monod(substrate, half_sat_const, exponent): ''' return (substrate / (half_sat_const + substrate)) ** exponent -# Temperature model (Arrhenius) +# Temperature model (Arrhenius) def temperature(temp, arr_a, arr_e): ''' - :param temp: temperature (will be imported from input excel file) [K] + :param temp: temperature (will be imported from input excel file) [K] :param arr_a: arrhenius constant (A) (Goldman et al., 1974) [unitless] :param arr_e: arrhenius exponential constant (E/R) (Goldman et al., 1974) [K] :return: temperature component of overall growth equation [unitless] ''' return arr_a * np.exp(-arr_e / temp) # Used equation from Goldman et al., 1974 -# Photoadaptation (_p1) +# Photoadaptation (_p1) def photoadaptation(i_avg, X_CHL, X_carbon, I_n, k_gamma): ''' :param i_avg: I, depth-averaged irradiance (calculated from 'attenuation' method) [uE/m^2/s] :param X_CHL: chlorophyll content of cells [g Chl/m^3] :param X_carbon: carbon content of cells (X_ALG + X_CH + X_LI) * i_C [g C/m^3] :param I_n: maximum incident PAR irradiance (“irradiance at noon”) [uE/m^2/s] - :param k_gamma: photoadaptation coefficient [unitless] - :return: photoadaptation rate [g Chl/m^3/d] + :param k_gamma: photoadaptation coefficient [unitless] + :return: photoadaptation rate [g Chl/m^3/d] ''' if X_carbon > 0: return 24 * ((0.2 * i_avg / I_n) / (k_gamma + (i_avg / I_n))) *\ @@ -266,19 +266,19 @@ def photoadaptation(i_avg, X_CHL, X_carbon, I_n, k_gamma): else: return 0 # Nutrients uptake (_p2, _p3, _p4, _p5, _p6) -def nutrient_uptake(X_ALG, quota, substrate, uptake_rate, half_sat_const, maximum_quota, subsistence_quota): +def nutrient_uptake(X_ALG, quota, substrate, uptake_rate, half_sat_const, maximum_quota, subsistence_quota): ''' :param X_ALG: algae biomass concentration (i.e., no storage products) [g COD/m^3] - :param quota: Q_N or Q_P [g N or g P/g COD] + :param quota: Q_N or Q_P [g N or g P/g COD] :param substrate: S_NH, S_NO or S_P [g N or g P/m^3] :param uptake_rate: V_NH, V_NO or V_P [g N or g P/g COD/d] :param half_sat_const: K_N or K_P [g N or g P/m^3] :param maximum_quota: Q_N_max or Q_P_max [g N or g P/g COD] - :param subsistence_quota: Q_N_min or Q_P_min [g N or g P/g COD] + :param subsistence_quota: Q_N_min or Q_P_min [g N or g P/g COD] :return: nutrient uptake rate [g N or g P/m^3/d] ''' return uptake_rate * monod(substrate, half_sat_const, 1) * X_ALG * \ - ((maximum_quota - quota) / (maximum_quota - subsistence_quota)) ** 0.01 + ((maximum_quota - quota) / (maximum_quota - subsistence_quota)) ** 0.01 # Maximum total photoautotrophic or heterotrophic-acetate or heterotrophic-glucose growth rate (_p7, _p10, _p11, _p15, _p18, _p19, _p23, _p26, _p27) def max_total_growth(X_ALG, mu_max, f_np, f_temp): @@ -287,7 +287,7 @@ def max_total_growth(X_ALG, mu_max, f_np, f_temp): :param mu_max: maximum specific growth rate [d^(-1)] :param f_np: inhibition factor by nitrogen or phosphorus (between 0 and 1) [unitless] :param f_temp: temperature correction factor (between 0 and 1) [unitless] - :return: maximum total growth rate for a particular mechanism, + :return: maximum total growth rate for a particular mechanism, without considering carbon source or light inhibition (= product of shared terms in growth-related equations) [g COD/m^3/d] ''' return mu_max * f_np * X_ALG * f_temp @@ -296,8 +296,8 @@ def max_total_growth(X_ALG, mu_max, f_np, f_temp): def growth_split(f_I, f_CH, f_LI, rho, Y_CH, Y_LI, K_STO): ''' :param f_I: irradiance response function (calculated from 'irrad_response' method) [unitless] - :param f_CH: ratio of stored carbohydrates to cells (X_CH / X_ALG) [g COD/g COD] - :param f_LI: ratio of stored lipids to cells (X_LI / X_ALG) [g COD/g COD] + :param f_CH: ratio of stored carbohydrates to cells (X_CH / X_ALG) [g COD/g COD] + :param f_LI: ratio of stored lipids to cells (X_LI / X_ALG) [g COD/g COD] :param rho: carbohydrate relative preference factor (calibrated in Guest et al., 2013) [unitless] :param Y_CH: yield of storage carbohydrates (as polyglucose, PG), Y_CH_PHO, Y_CH_NR_HET_ACE, or Y_CH_NR_HET_GLU [g COD/g COD] :param Y_LI: yield of storage lipids (as triacylglycerol, TAG), Y_LI_PHO, Y_LI_NR_HET_ACE, or Y_LI_NR_HET_GLU [g COD/g COD] @@ -308,7 +308,7 @@ def growth_split(f_I, f_CH, f_LI, rho, Y_CH, Y_LI, K_STO): numerators = np.asarray([K_STO * (1 - f_I), rho * f_CH, f_LI * Y_CH / Y_LI]) return numerators/(sum(numerators)) -# Part of storage equations (_p8, _p9, _p16, _p17, _p24, _p25) +# Part of storage equations (_p8, _p9, _p16, _p17, _p24, _p25) def storage_saturation(f, f_max, beta): ''' :param f: f_CH or f_LI [g COD/g COD] @@ -323,7 +323,7 @@ def max_total_maintenance(X_ALG, m_ATP): ''' :param X_ALG: algae biomass concentration (i.e., no storage products) [g COD/m^3] :param m_ATP: specific maintenance rate [g ATP/g COD/d] - :return: maximum total maintenance rate for a particular mechanism + :return: maximum total maintenance rate for a particular mechanism (= product of shared terms in maintenance-related equations) [g COD/m^3/d] ''' return m_ATP * X_ALG @@ -331,8 +331,8 @@ def max_total_maintenance(X_ALG, m_ATP): # Split the total maintenance rate between three processes (_p12, _p13, _p14, _p20, _p21, _p22, _p28, _p29, _p30) def maintenance_split(f_CH, f_LI, rho, Y_CH, Y_LI, Y_X_ALG, Y_ATP, K_STO): ''' - :param f_CH: ratio of stored carbohydrates to cells (X_CH / X_ALG) [g COD/g COD] - :param f_LI: ratio of stored lipids to cells (X_LI / X_ALG) [g COD/g COD] + :param f_CH: ratio of stored carbohydrates to cells (X_CH / X_ALG) [g COD/g COD] + :param f_LI: ratio of stored lipids to cells (X_LI / X_ALG) [g COD/g COD] :param rho: carbohydrate relative preference factor (calibrated in Guest et al., 2013) [unitless] :param Y_CH: yield of storage carbohydrates (as polyglucose, PG), Y_CH_PHO, Y_CH_NR_HET_ACE, or Y_CH_NR_HET_GLU [g COD/g COD] :param Y_LI: yield of storage lipids (as triacylglycerol, TAG), Y_LI_PHO, Y_LI_NR_HET_ACE, or Y_LI_NR_HET_GLU [g COD/g COD] @@ -358,8 +358,36 @@ def storage(X_ALG, f_np, response, saturation, storage_rate): ''' return storage_rate * saturation * (1 - f_np) * response * X_ALG -def rhos_pm2_asm2d(state_arr, params): - +'''added from asm2d''' + +# Hydrolysis (_p31, _p32, _p33) +def hydrolysis(X_S, X_H, K_h, K_X): + ''' + :param X_S: concentration of slowly biodegradable substrates + :param X_H: concentration of heterotrophic organisms (including denitrifer) + :param K_h: hydrolysis rate constant + :param K_X: slowly biodegradable substrate half saturation coefficient for hydrolysis + :return: shared parts of hydrolysis equations + ''' + return K_h * (X_S/X_H) / (K_X + X_S/X_H) * X_H + +# Growth in ASM2d (_p34, _p35, _p36, _p37, _p40) +def growth_asm2d(S_NH, S_P, S_ALK, mu, X, K_NH4, K_P, K_ALK): + ''' + :param S_NH: concentration of dissolved ammonium + :param S_P: concentration of dissolved phosphorus + :param S_ALK: concentration of alkalinity + :param mu: maximum specific growth rate (mu_H or mu_AUT) + :param X: concentration of biomass (X_H or X_AUT) + :param K_NH4: ammonium (nutrient) half saturation coefficient (K_NH4_H or K_NH4_AUT) + :param K_P: phosphorus (nutrient) half saturation coefficient (K_P_H or K_P_AUT) + :param K_ALK: alkalinity half saturation coefficient (K_ALK_H or K_ALK_AUT) + :return: shared parts of growth-related equations + ''' + return mu * S_NH/(K_NH4+S_NH) * S_P/(K_P+S_P) * S_ALK/(K_ALK + S_ALK) * X + +def rhos_pm2asm2d(state_arr, params): + # extract values of state variables c_arr = state_arr[:21] temp = state_arr[22] @@ -368,26 +396,23 @@ def rhos_pm2_asm2d(state_arr, params): # Q = state_arr[14] # Flow rate # t = state_arr[15] # time - # light = calc_irrad(t) # when to use calculated light (I_0) - - # X_CHL, X_ALG, X_CH, X_LI, S_CO2, S_A, S_F, S_O2, S_NH, S_NO, S_P, X_N_ALG, X_P_ALG, H2O = c_arr X_CHL, X_ALG, X_CH, X_LI, S_CO2, S_A, S_F, S_O2, S_NH, S_NO, S_P, X_N_ALG, X_P_ALG, S_N2, S_ALK, S_I, X_I, X_S, X_H, X_AUT, H2O = c_arr - + # extract values of parameters - cmps = params['cmps'] + cmps = params['cmps'] a_c = params['a_c'] I_n = params['I_n'] arr_a = params['arr_a'] arr_e = params['arr_e'] beta_1 = params['beta_1'] - beta_2 = params['beta_2'] + beta_2 = params['beta_2'] b_reactor = params['b_reactor'] I_opt = params['I_opt'] k_gamma = params['k_gamma'] K_N = params['K_N'] K_P = params['K_P'] K_A = params['K_A'] - K_F = params['K_F'] + K_F = params['K_F'] rho = params['rho'] K_STO = params['K_STO'] f_CH_max = params['f_CH_max'] @@ -417,7 +442,7 @@ def rhos_pm2_asm2d(state_arr, params): Y_LI_NR_HET_GLU = params['Y_LI_NR_HET_GLU'] Y_X_ALG_HET_GLU = params['Y_X_ALG_HET_GLU'] n_dark = params['n_dark'] - + '''added from asm2d''' f_SI = params['f_SI'] Y_H = params['Y_H'] @@ -452,17 +477,17 @@ def rhos_pm2_asm2d(state_arr, params): # intermediate variables f_CH = ratio(X_CH, X_ALG, 0, f_CH_max) f_LI = ratio(X_LI, X_ALG, 0, f_LI_max) - + # Q_N = ratio(X_N_ALG, X_ALG, Q_N_min, Q_N_max) # Q_P = ratio(X_P_ALG, X_ALG, Q_P_min, Q_P_max) - + alg_iN, alg_iP = cmps.X_ALG.i_N, cmps.X_ALG.i_P Q_N = ratio(X_N_ALG+X_ALG*alg_iN, X_ALG, Q_N_min, Q_N_max) Q_P = ratio(X_P_ALG+X_ALG*alg_iP, X_ALG, Q_P_min, Q_P_max) - + idx = cmps.indices(['X_ALG', 'X_CH', 'X_LI']) X_bio = np.array([X_ALG, X_CH, X_LI]) - X_TSS = sum(X_bio * cmps.i_mass[idx]) + X_TSS = sum(X_bio * cmps.i_mass[idx]) X_carbon = sum(X_bio * cmps.i_C[idx]) i_avg = attenuation(light, X_TSS, a_c, b_reactor) @@ -480,33 +505,33 @@ def rhos_pm2_asm2d(state_arr, params): max_total_growth_rho = max_total_growth(X_ALG, mu_max, f_np, f_temp) max_maintenance_rho = max_total_maintenance(X_ALG, m_ATP) # light = calc_irrad(t) - + # calculate kinetic rate values - rhos = np.empty(30) - + rhos = np.empty(41) + rhos[0] = photoadaptation(i_avg, X_CHL, X_carbon, I_n, k_gamma) - + rhos[1] = nutrient_uptake(X_ALG, Q_N, S_NH, V_NH, K_N, Q_N_max, Q_N_min) rhos[[2,3,4]] = nutrient_uptake(X_ALG, Q_N, S_NO, V_NO, K_N, Q_N_max, Q_N_min) * (K_N/(K_N + S_NH)) rhos[5] = nutrient_uptake(X_ALG, Q_P, S_P, V_P, K_P, Q_P_max, Q_P_min) - + rhos[[6,9,10]] = max_total_growth_rho \ - * growth_split(f_I, f_CH, f_LI, rho, Y_CH_PHO, Y_LI_PHO, K_STO) + * growth_split(f_I, f_CH, f_LI, rho, Y_CH_PHO, Y_LI_PHO, K_STO) rhos[6] *= f_I - rhos[[9,10]] *= dark_response - + rhos[[9,10]] *= dark_response + rhos[[14,17,18]] = max_total_growth_rho \ * acetate_response \ * growth_split(f_I, f_CH, f_LI, rho, Y_CH_NR_HET_ACE, Y_LI_NR_HET_ACE, K_STO) - + rhos[[22,25,26]] = max_total_growth_rho \ * glucose_response \ * growth_split(f_I, f_CH, f_LI, rho, Y_CH_NR_HET_GLU, Y_LI_NR_HET_GLU, K_STO) - + rhos[[11,12,13]] = max_maintenance_rho \ * maintenance_split(f_CH, f_LI, rho, Y_CH_PHO, Y_LI_PHO, Y_X_ALG_PHO, Y_ATP_PHO, K_STO) - + rhos[[19,20,21]] = max_maintenance_rho \ * maintenance_split(f_CH, f_LI, rho, Y_CH_NR_HET_ACE, Y_LI_NR_HET_ACE, Y_X_ALG_HET_ACE, Y_ATP_HET_ACE, K_STO) @@ -520,11 +545,25 @@ def rhos_pm2_asm2d(state_arr, params): rhos[23] = storage(X_ALG, f_np, glucose_response, f_sat_CH, q_CH) rhos[24] = storage(X_ALG, f_np, glucose_response, f_sat_LI, q_LI) * (f_CH / f_CH_max) + rhos[30] = hydrolysis(X_S, X_H, K_h, K_X) * monod(S_O2, K_O2, 1) + rhos[31] = hydrolysis(X_S, X_H, K_h, K_X) * eta_NO3 * monod(K_O2, S_O2, 1) * monod(S_NO, K_NO3, 1) + rhos[32] = hydrolysis(X_S, X_H, K_h, K_X) * eta_fe * monod(K_O2, S_O2, 1) * monod(K_NO3, S_NO, 1) + + rhos[33] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * monod(S_O2, K_O2_H, 1) * monod(S_F, K_F, 1) * monod(S_F, S_A, 1) + rhos[34] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * monod(S_O2, K_O2_H, 1) * monod(S_A, K_A_H, 1) * monod(S_A, S_F, 1) + rhos[35] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * eta_NO3_H * monod(K_O2_H, S_O2, 1) * monod(S_NO, K_NO3_H, 1) * monod(S_F, K_F, 1) * monod(S_F, S_A, 1) + rhos[36] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * eta_NO3_H * monod(K_O2_H, S_O2, 1) * monod(S_NO, K_NO3_H, 1) * monod(S_A, K_A_H, 1) * monod(S_A, S_F, 1) + rhos[37] = q_fe * monod(K_O2_H, S_O2, 1) * monod(K_NO3_H, S_NO, 1) * monod(S_F, K_fe, 1) * monod(S_ALK, K_ALK_H, 1) * X_H + rhos[38] = b_H * X_H + + rhos[39] = growth_asm2d(S_NH, S_P, S_ALK, mu_AUT, X_AUT, K_NH4_AUT, K_P_AUT, K_ALK_AUT) * monod(S_O2, K_O2_AUT, 1) + rhos[40] = b_AUT * X_AUT + return rhos #%% # ============================================================================= -# PM2 class +# PM2ASM2d class # ============================================================================= class PM2ASM2d(CompiledProcesses): @@ -533,82 +572,82 @@ class PM2ASM2d(CompiledProcesses): ---------- components: class:`CompiledComponents`, optional Components corresponding to each entry in the stoichiometry array, - defaults to thermosteam.settings.chemicals. + defaults to thermosteam.settings.chemicals. a_c : float, optional PAR absorption coefficient on a TSS (total suspended solids) basis, in [m^2/g TSS]. - The default is 0.049. + The default is 0.049. I_n : float, optional Maximum incident PAR irradiance (“irradiance at noon”), in [uE/m^2/s]. - The default is 250. + The default is 250. arr_a : float, optional Arrhenius constant (A), in [unitless]. - The default is 1.8 * 10**10. + The default is 1.8 * 10**10. arr_e : float, optional Arrhenius exponential constant (E/R), in [K]. - The default is 6842. + The default is 6842. beta_1 : float, optional Power coefficient for carbohydrate storage inhibition, in [unitless]. - The default is 2.90. + The default is 2.90. beta_2 : float, optional Power coefficient for lipid storage inhibition, in [unitless]. - The default is 3.50. + The default is 3.50. b_reactor : float, optional Thickness of reactor along light path, in [m]. - The default is 0.03. + The default is 0.03. I_opt : float, optional Optimal irradiance, in [uE/m^2/s]. - The default is 300. + The default is 300. k_gamma : float, optional Photoadaptation coefficient, in [unitless]. - The default is 0.00001. + The default is 0.00001. K_N : float, optional Nitrogen half-saturation constant, in [g N/m^3]. - The default is 0.1. + The default is 0.1. K_P : float, optional Phosphorus half-saturation constant, in [g P/m^3]. - The default is 1.0. + The default is 1.0. K_A : float, optional Organic carbon half-saturation constant (acetate) (Wagner, 2016), in [g COD/m^3]. - The default is 6.3. + The default is 6.3. K_F : float, optional Organic carbon half-saturation constant (glucose); assumes K_A = K_F, in [g COD/m^3]. - The default is 6.3. + The default is 6.3. rho : float, optional Carbohydrate relative preference factor (calibrated in Guest et al., 2013), in [unitless]. - The default is 1.186. + The default is 1.186. K_STO : float, optional Half-saturation constant for stored organic carbon (calibrated in Guest et al., 2013), in [g COD/g COD]. - The default is 1.566. + The default is 1.566. f_CH_max : float, optional Maximum achievable ratio of stored carbohydrates to functional cells, in [g COD/g COD]. - The default is 0.819. + The default is 0.819. f_LI_max : float, optional Maximum achievable ratio of stored lipids to functional cells, in [g COD/g COD]. - The default is 3.249. + The default is 3.249. m_ATP : float, optional Specific maintenance rate, in [g ATP/g COD/d]. - The default is 15.835. + The default is 15.835. mu_max : float, optional Maximum specific growth rate, in [d^(-1)]. - The default is 1.969. + The default is 1.969. q_CH : float, optional Maximum specific carbohydrate storage rate, in [g COD/g COD/d]. - The default is 0.594. + The default is 0.594. q_LI : float, optional Maximum specific lipid storage rate, in [g COD/g COD/d]. - The default is 0.910. + The default is 0.910. Q_N_max : float, optional Maximum nitrogen quota, in [g N/g COD]. - The default is 0.417. + The default is 0.417. Q_N_min : float, optional Nitrogen subsistence quota, in [g N/g COD]. - The default is 0.082. + The default is 0.082. Q_P_max : float, optional Maximum phosphorus quota, in [g P/g COD]. The default is 0.092. Q_P_min : float, optional Phosphorus subsistence quota; assumes N:P ratio of 5:1, in [g P/g COD]. - The default is 0.0163. + The default is 0.0163. V_NH : float, optional Maximum specific ammonium uptake rate (calibrated in Guest et al., 2013), in [g N/g COD/d]. The default is 0.254. @@ -626,10 +665,10 @@ class PM2ASM2d(CompiledProcesses): The default is 55.073. Y_CH_PHO : float, optional Yield of storage carbohydrate (as polyglucose, PG) on CO2 fixed to G3P, in [g COD/g CO2]. - The default is 0.754. + The default is 0.754. Y_LI_PHO : float, optional Yield of storage lipids (as triacylglycerol, TAG) on CO2 fixed to G3P, in [g COD/g CO2]. - The default is 0.901. + The default is 0.901. Y_X_ALG_PHO : float, optional Yield of carbon-accumulating phototrophic organisms on CO2 fixed to G3P, in [g COD/g CO2]. The default is 0.450. @@ -672,46 +711,89 @@ class PM2ASM2d(CompiledProcesses): n_dark: float, optional Dark growth reduction factor, in [unitless] The default is 0.7. + + + + + f_SI = params['f_SI'] + Y_H = params['Y_H'] + f_XI_H = params['f_XI_H'] + Y_A = params['Y_A'] + f_XI_AUT = params['f_XI_AUT'] + K_h = params['K_h'] + eta_NO3 = params['eta_NO3'] + eta_fe = params['eta_fe'] + K_O2 = params['K_O2'] + K_NO3 = params['K_NO3'] + K_X = params['K_X'] + mu_H = params['mu_H'] + q_fe = params['q_fe'] + eta_NO3_H = params['eta_NO3_H'] + b_H = params['b_H'] + K_O2_H = params['K_O2_H'] + K_F_H = params['K_F_H'] # K_F overlaps with PM2 -> change into K_F_H + K_fe = params['K_fe'] + K_A_H = params['K_A_H'] + K_NO3_H = params['K_NO3_H'] + K_NH4_H = params['K_NH4_H'] + K_P_H = params['K_P_H'] + K_ALK_H = params['K_ALK_H'] + mu_AUT = params['mu_AUT'] + b_AUT = params['b_AUT'] + K_O2_AUT = params['K_O2_AUT'] + K_NH4_AUT = params['K_NH4_AUT'] + K_ALK_AUT = params['K_ALK_AUT'] + K_P_AUT = params['K_P_AUT'] + + + + + + + path : str, optional - Alternative file path for the Petersen matrix. + Alternative file path for the Petersen matrix. The default is None. Examples -------- >>> from qsdsan import processes as pc - >>> cmps = pc.create_pm2_cmps() - >>> pm2 = pc.PM2() - >>> pm2.show() - PM2([photoadaptation, ammonium_uptake, nitrate_uptake_pho, nitrate_uptake_ace, nitrate_uptake_glu, phosphorus_uptake, - growth_pho, carbohydrate_storage_pho, lipid_storage_pho, carbohydrate_growth_pho, lipid_growth_pho, - carbohydrate_maintenance_pho, lipid_maintenance_pho, endogenous_respiration_pho, - growth_ace, carbohydrate_storage_ace, lipid_storage_ace, carbohydrate_growth_ace, lipid_growth_ace, - carbohydrate_maintenance_ace, lipid_maintenance_ace, endogenous_respiration_ace, - growth_glu, carbohydrate_storage_glu, lipid_storage_glu, carbohydrate_growth_glu, lipid_growth_glu, - carbohydrate_maintenance_glu, lipid_maintenance_glu, endogenous_respiration_glu]) + >>> cmps = pc.create_pm2asm2d_cmps() + >>> pm2asm2d = pc.PM2ASM2d() + >>> pm2asm2d.show() + PM2ASM2d([photoadaptation, ammonium_uptake, nitrate_uptake_pho, nitrate_uptake_ace, nitrate_uptake_glu, phosphorus_uptake, + growth_pho, carbohydrate_storage_pho, lipid_storage_pho, carbohydrate_growth_pho, lipid_growth_pho, + carbohydrate_maintenance_pho, lipid_maintenance_pho, endogenous_respiration_pho, + growth_ace, carbohydrate_storage_ace, lipid_storage_ace, carbohydrate_growth_ace, lipid_growth_ace, + carbohydrate_maintenance_ace, lipid_maintenance_ace, endogenous_respiration_ace, + growth_glu, carbohydrate_storage_glu, lipid_storage_glu, carbohydrate_growth_glu, lipid_growth_glu, + carbohydrate_maintenance_glu, lipid_maintenance_glu, endogenous_respiration_glu, + aero_hydrolysis, anox_hydrolysis, anae_hydrolysis, + hetero_growth_S_F, hetero_growth_S_A, denitri_S_F, denitri_S_A, ferment, hetero_lysis, + auto_aero_growth, auto_lysis]) ''' _shared_params = ('Y_CH_PHO', 'Y_LI_PHO', 'Y_X_ALG_PHO', 'Y_CH_NR_HET_ACE', 'Y_LI_NR_HET_ACE', 'Y_X_ALG_HET_ACE', 'Y_CH_NR_HET_GLU', 'Y_LI_NR_HET_GLU', 'Y_X_ALG_HET_GLU') - + _stoichio_params = ('Y_CH_ND_HET_ACE', 'Y_LI_ND_HET_ACE', 'Y_CH_ND_HET_GLU', 'Y_LI_ND_HET_GLU', *_shared_params) - - _kinetic_params = ('a_c', 'I_n', 'arr_a', 'arr_e', 'beta_1', 'beta_2', 'b_reactor', 'I_opt', 'k_gamma', - 'K_N', 'K_P', 'K_A', 'K_F', 'rho', 'K_STO', 'f_CH_max', 'f_LI_max', 'm_ATP', 'mu_max', - 'q_CH', 'q_LI', 'Q_N_max', 'Q_N_min', 'Q_P_max', 'Q_P_min', 'V_NH', 'V_NO', 'V_P', 'exponent', + + _kinetic_params = ('a_c', 'I_n', 'arr_a', 'arr_e', 'beta_1', 'beta_2', 'b_reactor', 'I_opt', 'k_gamma', + 'K_N', 'K_P', 'K_A', 'K_F', 'rho', 'K_STO', 'f_CH_max', 'f_LI_max', 'm_ATP', 'mu_max', + 'q_CH', 'q_LI', 'Q_N_max', 'Q_N_min', 'Q_P_max', 'Q_P_min', 'V_NH', 'V_NO', 'V_P', 'exponent', 'Y_ATP_PHO', 'Y_ATP_HET_ACE', 'Y_ATP_HET_GLU', *_shared_params, 'n_dark', 'cmps') def __new__(cls, components=None, - a_c=0.049, I_n=250, arr_a=1.8e10, arr_e=6842, beta_1=2.90, beta_2=3.50, b_reactor=0.03, I_opt=300, k_gamma=1e-5, - K_N=0.1, K_P=1.0, K_A=6.3, K_F=6.3, rho=1.186, K_STO=1.566, - f_CH_max=0.819, f_LI_max=3.249, m_ATP=15.835, mu_max=1.969, q_CH=0.594, q_LI=0.910, + a_c=0.049, I_n=250, arr_a=1.8e10, arr_e=6842, beta_1=2.90, beta_2=3.50, b_reactor=0.03, I_opt=300, k_gamma=1e-5, + K_N=0.1, K_P=1.0, K_A=6.3, K_F=6.3, rho=1.186, K_STO=1.566, + f_CH_max=0.819, f_LI_max=3.249, m_ATP=15.835, mu_max=1.969, q_CH=0.594, q_LI=0.910, Q_N_max=0.417, Q_N_min=0.082, Q_P_max=0.092, Q_P_min=0.0163, V_NH=0.254, V_NO=0.254, V_P=0.016, exponent=4, Y_ATP_PHO=55.073, Y_CH_PHO=0.754, Y_LI_PHO=0.901, Y_X_ALG_PHO=0.450, - Y_ATP_HET_ACE=39.623, Y_CH_NR_HET_ACE=0.625, Y_CH_ND_HET_ACE=0.600, + Y_ATP_HET_ACE=39.623, Y_CH_NR_HET_ACE=0.625, Y_CH_ND_HET_ACE=0.600, Y_LI_NR_HET_ACE=1.105, Y_LI_ND_HET_ACE=0.713, Y_X_ALG_HET_ACE=0.216, - Y_ATP_HET_GLU=58.114, Y_CH_NR_HET_GLU=0.917, Y_CH_ND_HET_GLU=0.880, + Y_ATP_HET_GLU=58.114, Y_CH_NR_HET_GLU=0.917, Y_CH_ND_HET_GLU=0.880, Y_LI_NR_HET_GLU=1.620, Y_LI_ND_HET_GLU=1.046, Y_X_ALG_HET_GLU=0.317, n_dark=0.7, path=None, **kwargs): @@ -721,7 +803,7 @@ def __new__(cls, components=None, conserved_for=('COD', 'C', 'N', 'P'), parameters=cls._stoichio_params, compile=False) - + if path == _path: _p3 = Process('nitrate_uptake_pho', 'S_NO -> [?]S_O2 + X_N_ALG', @@ -740,34 +822,34 @@ def __new__(cls, components=None, components=components, ref_component='X_N_ALG', conserved_for=('COD', 'C')) - + self.insert(2, _p3) self.insert(3, _p4) self.insert(4, _p5) self.compile(to_class=cls) - - self.set_rate_function(rhos_pm2_asm2d) - shared_values = (Y_CH_PHO, Y_LI_PHO, Y_X_ALG_PHO, - Y_CH_NR_HET_ACE, Y_LI_NR_HET_ACE, Y_X_ALG_HET_ACE, + + self.set_rate_function(rhos_pm2asm2d) + shared_values = (Y_CH_PHO, Y_LI_PHO, Y_X_ALG_PHO, + Y_CH_NR_HET_ACE, Y_LI_NR_HET_ACE, Y_X_ALG_HET_ACE, Y_CH_NR_HET_GLU, Y_LI_NR_HET_GLU, Y_X_ALG_HET_GLU) stoichio_values = (Y_CH_ND_HET_ACE, Y_LI_ND_HET_ACE, Y_CH_ND_HET_GLU, Y_LI_ND_HET_GLU, *shared_values) Q_N_min = max(self.Th_Q_N_min, Q_N_min) Q_P_min = max(self.Th_Q_P_min, Q_P_min) - kinetic_values = (a_c, I_n, arr_a, arr_e, beta_1, beta_2, b_reactor, I_opt, k_gamma, - K_N, K_P, K_A, K_F, rho, K_STO, f_CH_max, f_LI_max, m_ATP, mu_max, - q_CH, q_LI, Q_N_max, Q_N_min, Q_P_max, Q_P_min, V_NH, V_NO, V_P, exponent, + kinetic_values = (a_c, I_n, arr_a, arr_e, beta_1, beta_2, b_reactor, I_opt, k_gamma, + K_N, K_P, K_A, K_F, rho, K_STO, f_CH_max, f_LI_max, m_ATP, mu_max, + q_CH, q_LI, Q_N_max, Q_N_min, Q_P_max, Q_P_min, V_NH, V_NO, V_P, exponent, Y_ATP_PHO, Y_ATP_HET_ACE, Y_ATP_HET_GLU, *shared_values, n_dark, self._components) - + dct = self.__dict__ dct.update(kwargs) dct['_parameters'] = dict(zip(cls._stoichio_params, stoichio_values)) self.rate_function._params = dict(zip(cls._kinetic_params, kinetic_values)) return self - + def set_parameters(self, **parameters): '''Set values to stoichiometric and/or kinetic parameters.''' stoichio_only = {k:v for k,v in parameters.items() if k in self._stoichio_params} @@ -775,11 +857,11 @@ def set_parameters(self, **parameters): if self._stoichio_lambdified is not None: self.__dict__['_stoichio_lambdified'] = None if 'Q_N_min' in parameters.keys(): - if parameters['Q_N_min'] < self.Th_Q_N_min: + if parameters['Q_N_min'] < self.Th_Q_N_min: raise ValueError(f'Value for Q_N_min must not be less than the ' f'theoretical minimum {self.Th_Q_N_min}') if 'Q_P_min' in parameters.keys(): - if parameters['Q_P_min'] < self.Th_Q_P_min: + if parameters['Q_P_min'] < self.Th_Q_P_min: raise ValueError(f'Value for Q_P_min must not be less than the ' f'theoretical minimum {self.Th_Q_P_min}') self.rate_function.set_param(**parameters) @@ -787,7 +869,7 @@ def set_parameters(self, **parameters): @property def Th_Q_N_min(self): return abs(self.stoichiometry.loc['growth_pho', 'X_N_ALG'])*1.001 - + @property def Th_Q_P_min(self): return abs(self.stoichiometry.loc['growth_pho', 'X_P_ALG'])*1.001 \ No newline at end of file From 293e0a6d633124823ed2beb7e2a465c9b89dfd01 Mon Sep 17 00:00:00 2001 From: Ga-Yeong Kim <47093338+GaYeongKim@users.noreply.github.com> Date: Tue, 23 Apr 2024 15:25:23 -0500 Subject: [PATCH 03/32] Update _pm2asm2d.py --- qsdsan/processes/_pm2asm2d.py | 168 ++++++++++++++++++++++++---------- 1 file changed, 121 insertions(+), 47 deletions(-) diff --git a/qsdsan/processes/_pm2asm2d.py b/qsdsan/processes/_pm2asm2d.py index 7a0d98b2..7b5adee6 100644 --- a/qsdsan/processes/_pm2asm2d.py +++ b/qsdsan/processes/_pm2asm2d.py @@ -393,8 +393,8 @@ def rhos_pm2asm2d(state_arr, params): temp = state_arr[22] light = state_arr[23] # imported from input file assumed - # Q = state_arr[14] # Flow rate - # t = state_arr[15] # time + # Q = state_arr[21] # Flow rate + # t = state_arr[22] # time X_CHL, X_ALG, X_CH, X_LI, S_CO2, S_A, S_F, S_O2, S_NH, S_NO, S_P, X_N_ALG, X_P_ALG, S_N2, S_ALK, S_I, X_I, X_S, X_H, X_AUT, H2O = c_arr @@ -710,47 +710,94 @@ class PM2ASM2d(CompiledProcesses): The default is 0.317. n_dark: float, optional Dark growth reduction factor, in [unitless] - The default is 0.7. - - - - - f_SI = params['f_SI'] - Y_H = params['Y_H'] - f_XI_H = params['f_XI_H'] - Y_A = params['Y_A'] - f_XI_AUT = params['f_XI_AUT'] - K_h = params['K_h'] - eta_NO3 = params['eta_NO3'] - eta_fe = params['eta_fe'] - K_O2 = params['K_O2'] - K_NO3 = params['K_NO3'] - K_X = params['K_X'] - mu_H = params['mu_H'] - q_fe = params['q_fe'] - eta_NO3_H = params['eta_NO3_H'] - b_H = params['b_H'] - K_O2_H = params['K_O2_H'] - K_F_H = params['K_F_H'] # K_F overlaps with PM2 -> change into K_F_H - K_fe = params['K_fe'] - K_A_H = params['K_A_H'] - K_NO3_H = params['K_NO3_H'] - K_NH4_H = params['K_NH4_H'] - K_P_H = params['K_P_H'] - K_ALK_H = params['K_ALK_H'] - mu_AUT = params['mu_AUT'] - b_AUT = params['b_AUT'] - K_O2_AUT = params['K_O2_AUT'] - K_NH4_AUT = params['K_NH4_AUT'] - K_ALK_AUT = params['K_ALK_AUT'] - K_P_AUT = params['K_P_AUT'] - - - - - - - + The default is 0.7. + f_SI : float, optional + Production of soluble inerts in hydrolysis, in [g COD/g COD]. + The default is 0.0. + Y_H : float, optional + Heterotrophic yield coefficient, in[g COD/g COD]. + The default is 0.625. + f_XI_H : float, optional + Fraction of inert COD generated in heterotrophic biomass lysis, in [g COD/g COD]. + The default is 0.1. + Y_A : float, optional + Autotrophic yield, in [g COD/g N]. + The default is 0.24. + f_XI_AUT : float, optional + Fraction of inert COD generated in autotrophic biomass lysis, in [g COD/g COD]. + The default is 0.1. + K_h : float, optional + Hydrolysis rate constant, in [d^(-1)]. + The default is 3.0. + eta_NO3 : float, optional + Reduction factor for anoxic hydrolysis, dimensionless. + The default is 0.6. + eta_fe : float, optional + Anaerobic hydrolysis reduction factor, dimensionless. + The default is 0.4. + K_O2 : float, optional + O2 half saturation coefficient for hydrolysis, in [g O2/m^3]. + The default is 0.2. + K_NO3 : float, optional + Nitrate half saturation coefficient for hydrolysis, in [g N/m^3]. + The default is 0.5. + K_X : float, optional + Slowly biodegradable substrate half saturation coefficient for hydrolysis, in [g COD/g COD]. + The default is 0.1. + mu_H : float, optional + Heterotrophic maximum specific growth rate, in [d^(-1)]. + The default is 6.0. + q_fe : float, optional + Fermentation maximum rate, in [d^(-1)]. + The default is 3.0. + eta_NO3_H : float, optional + Reduction factor for anoxic heterotrophic growth, dimensionless. + The default is 0.8. + b_H : float, optional + Lysis and decay rate constant, in [d^(-1)]. + The default is 0.4. + K_O2_H : float, optional + O2 half saturation coefficient for heterotrophs, in [g O2/m^3]. + The default is 0.2. + K_F_H : float, optional + Fermentable substrate half saturation coefficient for heterotrophic growth (K_F in ASM2d), in [g COD/m^3]. + The default is 4.0. + K_fe : float, optional + Fermentable substrate half saturation coefficient for fermentation, in [g COD/m^3]. + The default is 4.0. + K_A_H : float, optional + VFA half saturation coefficient for heterotrophs, in [g COD/m^3]. + The default is 4.0. + K_NO3_H : float, optional + Nitrate half saturation coefficient for heterotrophs, in [g N/m^3]. + The default is 0.5. + K_NH4_H : float, optional + Ammonium (nutrient) half saturation coefficient for heterotrophs, in [g N/m^3]. + The default is 0.05. + K_P_H : float, optional + Phosphorus (nutrient) half saturation coefficient for heterotrophs, in [g P/m^3]. + The default is 0.01. + K_ALK_H : float, optional + Alkalinity half saturation coefficient for heterotrophs, in [mol(HCO3-)/m^3]. (user input unit, converted as C) + The default is 0.1. + mu_AUT : float, optional + Autotrophic maximum specific growth rate, in [d^(-1)]. + The default is 1.0. + b_AUT : float, optional + Autotrophic decay rate, in [d^(-1)]. + The default is 0.15. + K_O2_AUT : float, optional + O2 half saturation coefficient for autotrophs, in [g O2/m^3]. + The default is 0.5. + K_NH4_AUT : float, optional + Ammonium (nutrient) half saturation coefficient for autotrophs, in [g N/m^3]. + The default is 1.0. + K_ALK_AUT : float, optional + Alkalinity half saturation coefficient for autotrophs, in [mol(HCO3-)/m^3]. (user input unit, converted as C) + The default is 0.5. + K_P_AUT : float, optional + Phosphorus (nutrient) half saturation coefficient for autotrophs, in [g P/m^3]. + The default is 0.01. path : str, optional Alternative file path for the Petersen matrix. The default is None. @@ -771,6 +818,17 @@ class PM2ASM2d(CompiledProcesses): aero_hydrolysis, anox_hydrolysis, anae_hydrolysis, hetero_growth_S_F, hetero_growth_S_A, denitri_S_F, denitri_S_A, ferment, hetero_lysis, auto_aero_growth, auto_lysis]) + + References + ---------- + .. [1] Henze, M.; Gujer, W.; Mino, T.; Loosdrecht, M. van. Activated Sludge + Models: ASM1, ASM2, ASM2d and ASM3; IWA task group on mathematical modelling + for design and operation of biological wastewater treatment, Ed.; IWA + Publishing: London, 2000. + .. [2] Rieger, L.; Gillot, S.; Langergraber, G.; Ohtsuki, T.; Shaw, A.; Takács, + I.; Winkler, S. Guidelines for Using Activated Sludge Models; IWA Publishing: + London, New York, 2012; Vol. 11. + https://doi.org/10.2166/9781780401164. ''' _shared_params = ('Y_CH_PHO', 'Y_LI_PHO', 'Y_X_ALG_PHO', @@ -778,12 +836,16 @@ class PM2ASM2d(CompiledProcesses): 'Y_CH_NR_HET_GLU', 'Y_LI_NR_HET_GLU', 'Y_X_ALG_HET_GLU') _stoichio_params = ('Y_CH_ND_HET_ACE', 'Y_LI_ND_HET_ACE', 'Y_CH_ND_HET_GLU', 'Y_LI_ND_HET_GLU', + 'f_SI', 'Y_H', 'f_XI_H', 'Y_A', 'f_XI_AUT', *_shared_params) _kinetic_params = ('a_c', 'I_n', 'arr_a', 'arr_e', 'beta_1', 'beta_2', 'b_reactor', 'I_opt', 'k_gamma', 'K_N', 'K_P', 'K_A', 'K_F', 'rho', 'K_STO', 'f_CH_max', 'f_LI_max', 'm_ATP', 'mu_max', 'q_CH', 'q_LI', 'Q_N_max', 'Q_N_min', 'Q_P_max', 'Q_P_min', 'V_NH', 'V_NO', 'V_P', 'exponent', - 'Y_ATP_PHO', 'Y_ATP_HET_ACE', 'Y_ATP_HET_GLU', *_shared_params, 'n_dark', 'cmps') + 'Y_ATP_PHO', 'Y_ATP_HET_ACE', 'Y_ATP_HET_GLU', *_shared_params, 'n_dark', 'cmps', + 'K_h', 'eta_NO3', 'eta_fe', 'K_O2', 'K_NO3', 'K_X', 'mu_H', 'q_fe', 'eta_NO3_H', + 'b_H', 'K_O2_H', 'K_F_H', 'K_fe', 'K_A_H', 'K_NO3_H', 'K_NH4_H', 'K_P_H', + 'K_ALK_H', 'mu_AUT', 'b_AUT', 'K_O2_AUT', 'K_NH4_AUT', 'K_ALK_AUT', 'K_P_AUT') def __new__(cls, components=None, a_c=0.049, I_n=250, arr_a=1.8e10, arr_e=6842, beta_1=2.90, beta_2=3.50, b_reactor=0.03, I_opt=300, k_gamma=1e-5, @@ -795,12 +857,19 @@ def __new__(cls, components=None, Y_LI_NR_HET_ACE=1.105, Y_LI_ND_HET_ACE=0.713, Y_X_ALG_HET_ACE=0.216, Y_ATP_HET_GLU=58.114, Y_CH_NR_HET_GLU=0.917, Y_CH_ND_HET_GLU=0.880, Y_LI_NR_HET_GLU=1.620, Y_LI_ND_HET_GLU=1.046, Y_X_ALG_HET_GLU=0.317, n_dark=0.7, + f_SI=0.0, Y_H=0.625, f_XI_H=0.1, Y_A=0.24, f_XI_AUT=0.1, + K_h=3.0, eta_NO3=0.6, eta_fe=0.4, K_O2=0.2, K_NO3=0.5, K_X=0.1, + mu_H=6.0, q_fe=3.0, eta_NO3_H=0.8, b_H=0.4, K_O2_H=0.2, K_F_H=4.0, + K_fe=4.0, K_A_H=4.0, K_NO3_H=0.5, K_NH4_H=0.05, K_P_H=0.01, K_ALK_H=0.1, + mu_AUT=1.0, b_AUT=0.15, K_O2_AUT=0.5, K_NH4_AUT=1.0, K_ALK_AUT=0.5, K_P_AUT=0.01, path=None, **kwargs): if not path: path = _path + self = Processes.load_from_file(path, components=components, - conserved_for=('COD', 'C', 'N', 'P'), + conserved_for=('COD', 'C', 'N', 'P', 'charge'), + #conserved_for=('COD', 'C', 'N', 'P'), parameters=cls._stoichio_params, compile=False) @@ -834,14 +903,19 @@ def __new__(cls, components=None, Y_CH_NR_HET_ACE, Y_LI_NR_HET_ACE, Y_X_ALG_HET_ACE, Y_CH_NR_HET_GLU, Y_LI_NR_HET_GLU, Y_X_ALG_HET_GLU) stoichio_values = (Y_CH_ND_HET_ACE, Y_LI_ND_HET_ACE, Y_CH_ND_HET_GLU, Y_LI_ND_HET_GLU, - *shared_values) + f_SI, Y_H, f_XI_H, Y_A, f_XI_AUT, + *shared_values) Q_N_min = max(self.Th_Q_N_min, Q_N_min) Q_P_min = max(self.Th_Q_P_min, Q_P_min) kinetic_values = (a_c, I_n, arr_a, arr_e, beta_1, beta_2, b_reactor, I_opt, k_gamma, K_N, K_P, K_A, K_F, rho, K_STO, f_CH_max, f_LI_max, m_ATP, mu_max, q_CH, q_LI, Q_N_max, Q_N_min, Q_P_max, Q_P_min, V_NH, V_NO, V_P, exponent, Y_ATP_PHO, Y_ATP_HET_ACE, Y_ATP_HET_GLU, - *shared_values, n_dark, self._components) + *shared_values, n_dark, self._components, + K_h, eta_NO3, eta_fe, K_O2, K_NO3, K_X, mu_H, q_fe, eta_NO3_H, + b_H, K_O2_H, K_F_H, K_fe, K_A_H, K_NO3_H, K_NH4_H, K_P_H, + K_ALK_H*12, mu_AUT, b_AUT, K_O2_AUT, K_NH4_AUT, K_ALK_AUT*12, K_P_AUT, + ) dct = self.__dict__ dct.update(kwargs) From d39e31c21d45508e0a48a6a065e8aa786bd9e075 Mon Sep 17 00:00:00 2001 From: Ga-Yeong Kim <47093338+GaYeongKim@users.noreply.github.com> Date: Tue, 23 Apr 2024 17:02:10 -0500 Subject: [PATCH 04/32] split pm2 & asm2d stoichiometry path --- qsdsan/data/process_data/_pm2asm2d_1.tsv | 28 ++++++++++++++++++++++++ qsdsan/data/process_data/_pm2asm2d_2.tsv | 12 ++++++++++ qsdsan/processes/__init__.py | 4 ++++ qsdsan/processes/_pm2asm2d.py | 14 +++++++++--- 4 files changed, 55 insertions(+), 3 deletions(-) create mode 100644 qsdsan/data/process_data/_pm2asm2d_1.tsv create mode 100644 qsdsan/data/process_data/_pm2asm2d_2.tsv diff --git a/qsdsan/data/process_data/_pm2asm2d_1.tsv b/qsdsan/data/process_data/_pm2asm2d_1.tsv new file mode 100644 index 00000000..f3a977f5 --- /dev/null +++ b/qsdsan/data/process_data/_pm2asm2d_1.tsv @@ -0,0 +1,28 @@ + X_CHL X_ALG X_CH X_LI S_CO2 S_A S_F S_O2 S_NH S_NO S_P X_N_ALG X_P_ALG S_N2 S_ALK S_I X_I X_S X_H X_AUT +photoadaptation 1 +ammonium_uptake -1 1 +phosphorus_uptake -1 1 +growth_pho 1 ? 1 ? ? +carbohydrate_storage_pho 1 ? 1 +lipid_storage_pho 1 ? 1 +carbohydrate_growth_pho 1 (-Y_CH_PHO/Y_X_ALG_PHO) ? ? ? ? +lipid_growth_pho 1 (-Y_LI_PHO/Y_X_ALG_PHO) ? ? ? ? +carbohydrate_maintenance_pho -1 ? -1 +lipid_maintenance_pho -1 ? -1 +endogenous_respiration_pho -1 ? -1 ? ? +growth_ace 1 ? (-1)/Y_X_ALG_HET_ACE ? ? ? +carbohydrate_storage_ace 1 ? (-1)/Y_CH_ND_HET_ACE ? +lipid_storage_ace 1 ? (-1)/Y_LI_ND_HET_ACE ? +carbohydrate_growth_ace 1 (-Y_CH_NR_HET_ACE/Y_X_ALG_HET_ACE) ? ? ? ? +lipid_growth_ace 1 (-Y_LI_NR_HET_ACE/Y_X_ALG_HET_ACE) ? ? ? ? +carbohydrate_maintenance_ace -1 ? -1 +lipid_maintenance_ace -1 ? -1 +endogenous_respiration_ace -1 ? -1 ? ? +growth_glu 1 ? (-1)/Y_X_ALG_HET_GLU ? ? ? +carbohydrate_storage_glu 1 ? (-1)/Y_CH_ND_HET_GLU ? +lipid_storage_glu 1 ? (-1)/Y_LI_ND_HET_GLU ? +carbohydrate_growth_glu 1 (-Y_CH_NR_HET_GLU/Y_X_ALG_HET_GLU) ? ? ? ? +lipid_growth_glu 1 (-Y_LI_NR_HET_GLU/Y_X_ALG_HET_GLU) ? ? ? ? +carbohydrate_maintenance_glu -1 ? -1 +lipid_maintenance_glu -1 ? -1 +endogenous_respiration_glu -1 ? -1 ? ? diff --git a/qsdsan/data/process_data/_pm2asm2d_2.tsv b/qsdsan/data/process_data/_pm2asm2d_2.tsv new file mode 100644 index 00000000..8941b263 --- /dev/null +++ b/qsdsan/data/process_data/_pm2asm2d_2.tsv @@ -0,0 +1,12 @@ + X_CHL X_ALG X_CH X_LI S_CO2 S_A S_F S_O2 S_NH S_NO S_P X_N_ALG X_P_ALG S_N2 S_ALK S_I X_I X_S X_H X_AUT +aero_hydrolysis 1-f_SI ? ? ? f_SI -1 +anox_hydrolysis 1-f_SI ? ? ? f_SI -1 +anae_hydrolysis 1-f_SI ? ? ? f_SI -1 +hetero_growth_S_F (-1)/Y_H 1-1/Y_H ? ? ? 1 +hetero_growth_S_A (-1)/Y_H 1-1/Y_H ? ? ? 1 +denitri_S_F (-1)/Y_H ? (Y_H-1)/(20/7*Y_H) ? (1-Y_H)/(20/7*Y_H) ? 1 +denitri_S_A (-1)/Y_H ? (Y_H-1)/(20/7*Y_H) ? (1-Y_H)/(20/7*Y_H) ? 1 +ferment 1 -1 ? ? ? +hetero_lysis ? ? ? f_XI_H 1-f_XI_H -1 +auto_aero_growth (Y_A-32/7)/Y_A ? 1/Y_A ? ? 1 +auto_lysis ? ? ? f_XI_AUT 1-f_XI_AUT -1 diff --git a/qsdsan/processes/__init__.py b/qsdsan/processes/__init__.py index 9e0c5e9d..ddcf5d2c 100644 --- a/qsdsan/processes/__init__.py +++ b/qsdsan/processes/__init__.py @@ -21,6 +21,7 @@ from ._decay import * from ._kinetic_reaction import * from ._pm2 import * +from ._pm2asm2d import * from . import ( _aeration, @@ -30,6 +31,8 @@ _madm1, _decay, _kinetic_reaction, + _pm2, + _pm2asm2d, ) __all__ = ( @@ -41,4 +44,5 @@ *_decay.__all__, *_kinetic_reaction.__all__, *_pm2.__all__, + *_pm2asm2d.__all__, ) \ No newline at end of file diff --git a/qsdsan/processes/_pm2asm2d.py b/qsdsan/processes/_pm2asm2d.py index 7b5adee6..d56cdbd8 100644 --- a/qsdsan/processes/_pm2asm2d.py +++ b/qsdsan/processes/_pm2asm2d.py @@ -18,7 +18,9 @@ __all__ = ('create_pm2asm2d_cmps', 'PM2ASM2d') -_path = ospath.join(data_path, 'process_data/_pm2asm2d.tsv') +_path = ospath.join(data_path, 'process_data/_pm2asm2d_1.tsv') +_path_2 = ospath.join(data_path, 'process_data/_pm2asm2d_2.tsv') + # _load_components = settings.get_default_chemicals #%% @@ -868,10 +870,16 @@ def __new__(cls, components=None, self = Processes.load_from_file(path, components=components, - conserved_for=('COD', 'C', 'N', 'P', 'charge'), - #conserved_for=('COD', 'C', 'N', 'P'), + conserved_for=('COD', 'C', 'N', 'P'), + parameters=cls._stoichio_params, + compile=False) + + asm2d_processes = Processes.load_from_file(_path_2, + components=components, + conserved_for=('COD', 'N', 'P', 'charge'), parameters=cls._stoichio_params, compile=False) + self.extend(asm2d_processes) if path == _path: _p3 = Process('nitrate_uptake_pho', From 4ec854adadc3a6f0ed852f03137448078ff3a0ac Mon Sep 17 00:00:00 2001 From: Ga-Yeong Kim <47093338+GaYeongKim@users.noreply.github.com> Date: Sat, 27 Apr 2024 18:57:51 -0500 Subject: [PATCH 05/32] Update _pm2asm2d.py --- qsdsan/processes/_pm2asm2d.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/qsdsan/processes/_pm2asm2d.py b/qsdsan/processes/_pm2asm2d.py index d56cdbd8..223c10bd 100644 --- a/qsdsan/processes/_pm2asm2d.py +++ b/qsdsan/processes/_pm2asm2d.py @@ -446,11 +446,11 @@ def rhos_pm2asm2d(state_arr, params): n_dark = params['n_dark'] '''added from asm2d''' - f_SI = params['f_SI'] - Y_H = params['Y_H'] - f_XI_H = params['f_XI_H'] - Y_A = params['Y_A'] - f_XI_AUT = params['f_XI_AUT'] + # f_SI = params['f_SI'] + # Y_H = params['Y_H'] + # f_XI_H = params['f_XI_H'] + # Y_A = params['Y_A'] + # f_XI_AUT = params['f_XI_AUT'] K_h = params['K_h'] eta_NO3 = params['eta_NO3'] eta_fe = params['eta_fe'] @@ -551,9 +551,9 @@ def rhos_pm2asm2d(state_arr, params): rhos[31] = hydrolysis(X_S, X_H, K_h, K_X) * eta_NO3 * monod(K_O2, S_O2, 1) * monod(S_NO, K_NO3, 1) rhos[32] = hydrolysis(X_S, X_H, K_h, K_X) * eta_fe * monod(K_O2, S_O2, 1) * monod(K_NO3, S_NO, 1) - rhos[33] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * monod(S_O2, K_O2_H, 1) * monod(S_F, K_F, 1) * monod(S_F, S_A, 1) + rhos[33] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * monod(S_O2, K_O2_H, 1) * monod(S_F, K_F_H, 1) * monod(S_F, S_A, 1) rhos[34] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * monod(S_O2, K_O2_H, 1) * monod(S_A, K_A_H, 1) * monod(S_A, S_F, 1) - rhos[35] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * eta_NO3_H * monod(K_O2_H, S_O2, 1) * monod(S_NO, K_NO3_H, 1) * monod(S_F, K_F, 1) * monod(S_F, S_A, 1) + rhos[35] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * eta_NO3_H * monod(K_O2_H, S_O2, 1) * monod(S_NO, K_NO3_H, 1) * monod(S_F, K_F_H, 1) * monod(S_F, S_A, 1) rhos[36] = growth_asm2d(S_NH, S_P, S_ALK, mu_H, X_H, K_NH4_H, K_P_H, K_ALK_H) * eta_NO3_H * monod(K_O2_H, S_O2, 1) * monod(S_NO, K_NO3_H, 1) * monod(S_A, K_A_H, 1) * monod(S_A, S_F, 1) rhos[37] = q_fe * monod(K_O2_H, S_O2, 1) * monod(K_NO3_H, S_NO, 1) * monod(S_F, K_fe, 1) * monod(S_ALK, K_ALK_H, 1) * X_H rhos[38] = b_H * X_H From 9e4bd182a3f5bcb021cf7224c771d3b4ed102f0e Mon Sep 17 00:00:00 2001 From: Yalin Date: Fri, 12 Jul 2024 15:17:33 -0400 Subject: [PATCH 06/32] enable remote workflow --- .github/workflows/build-only.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/build-only.yml b/.github/workflows/build-only.yml index 5b5183e7..e8de2db5 100644 --- a/.github/workflows/build-only.yml +++ b/.github/workflows/build-only.yml @@ -4,9 +4,9 @@ name: build-only on: push: - branches: [beta, dev] + branches: [beta, dev, trial] pull_request: - branches: [beta, dev] + branches: [beta, dev, trial] jobs: build: From ceca243057082927b80cb84d01ad1176a9f62cb8 Mon Sep 17 00:00:00 2001 From: Yalin Date: Fri, 12 Jul 2024 15:19:01 -0400 Subject: [PATCH 07/32] try newer numpy --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 0a261e4c..c9e85842 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,7 +6,7 @@ SALib>=1.4.5 seaborn sympy>=1.8 matplotlib<=3.6.0 -numpy<2.0.0 +# numpy<2.0.0 # Specifically for docs sphinx From bcfdc1a180967a726aaed727d0e5fb494a307087 Mon Sep 17 00:00:00 2001 From: Yalin Date: Thu, 17 Oct 2024 12:43:24 -0400 Subject: [PATCH 08/32] more flexible coding on Reactor --- qsdsan/sanunits/_reactor.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/qsdsan/sanunits/_reactor.py b/qsdsan/sanunits/_reactor.py index d1da7fe7..89d5cec1 100644 --- a/qsdsan/sanunits/_reactor.py +++ b/qsdsan/sanunits/_reactor.py @@ -130,12 +130,13 @@ def _design(self): # if auxiliary, in our system, only K/O drum whose N, and V are provided # do not need to deal with self.F_vol_in (auxiliary unit has trouble doing this) ins_F_vol = self.F_vol_in + cmps = self.components + gas_IDs = [i.ID for i in cmps if i.locked_state=='g'] + sol_IDs = [i.ID for i in cmps if i.locked_state=='s'] for i in range(len(self.ins)): - ins_F_vol -= (self.ins[i].ivol['H2'] +\ - self.ins[i].ivol['CHG_catalyst'] +\ - self.ins[i].ivol['HT_catalyst'] +\ - self.ins[i].ivol['HC_catalyst']) - # not include gas (e.g. H2) + ins_F_vol -= sum(self.ins[i].ivol[gas_IDs]) + ins_F_vol -= sum(self.ins[i].ivol[sol_IDs]) + # not include gas (e.g. H2) and solids (e.g., catalysts) V_total = ins_F_vol * self.tau / self.V_wf P = self.P * _Pa_to_psi # Pa to psi length_to_diameter = self.length_to_diameter From beda37865c08308a63e5291e38542d0856623712 Mon Sep 17 00:00:00 2001 From: Yalin Date: Thu, 17 Oct 2024 20:44:40 -0400 Subject: [PATCH 09/32] fix bug related to `SanUnit.add_OPEX` --- qsdsan/_sanunit.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qsdsan/_sanunit.py b/qsdsan/_sanunit.py index 0cca88b1..23a5adfb 100644 --- a/qsdsan/_sanunit.py +++ b/qsdsan/_sanunit.py @@ -634,7 +634,7 @@ def results(self, with_units=True, include_utilities=True, results = super().results(with_units, include_utilities, include_total_cost, include_installed_cost, include_zeros, external_utilities, key_hook) - if not self.add_OPEX: self.add_OPEX = {'Additional OPEX': 0} + if not hasattr(self, 'add_OPEX'): self.add_OPEX = {'Additional OPEX': 0} for k, v in self.add_OPEX.items(): if not with_units: results.loc[(k, '')] = v From 838cc63634d1aca5bb88fbbe853e4eb295ac68d0 Mon Sep 17 00:00:00 2001 From: Yalin Date: Thu, 17 Oct 2024 20:44:59 -0400 Subject: [PATCH 10/32] add new inherited units from biosteam --- qsdsan/sanunits/_distillation.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/qsdsan/sanunits/_distillation.py b/qsdsan/sanunits/_distillation.py index b81cda4a..8f5040a6 100644 --- a/qsdsan/sanunits/_distillation.py +++ b/qsdsan/sanunits/_distillation.py @@ -20,6 +20,8 @@ __all__ = ( 'BinaryDistillation', 'ShortcutColumn', + 'MESHDistillation', + 'AdiabaticMultiStageVLEColumn', ) _lb_to_kg = qs.utils.auom('lb').conversion_factor('kg') @@ -53,6 +55,24 @@ class ShortcutColumn(bst.units.ShortcutColumn, qs.SanUnit): ''' biosteam.units.ShortcutColumn with QSDsan properties. + See Also + -------- + `biosteam.units.ShortcutColumn `_ + ''' + +class MESHDistillation(bst.units.MESHDistillation, qs.SanUnit): + ''' + biosteam.units.MESHDistillation with QSDsan properties. + + See Also + -------- + `biosteam.units.ShortcutColumn `_ + ''' + +class AdiabaticMultiStageVLEColumn(bst.units.AdiabaticMultiStageVLEColumn, qs.SanUnit): + ''' + biosteam.units.AdiabaticMultiStageVLEColumn with QSDsan properties. + See Also -------- `biosteam.units.ShortcutColumn `_ From ea1e7909932b1fb5247e7be22eee576056d2f7bf Mon Sep 17 00:00:00 2001 From: Yalin Date: Fri, 18 Oct 2024 10:28:31 -0400 Subject: [PATCH 11/32] more flexible setting in `Hydrocracking` --- qsdsan/sanunits/_hydroprocessing.py | 109 +++++++++++++++++++--------- 1 file changed, 73 insertions(+), 36 deletions(-) diff --git a/qsdsan/sanunits/_hydroprocessing.py b/qsdsan/sanunits/_hydroprocessing.py index d5e8601d..ce89cef4 100644 --- a/qsdsan/sanunits/_hydroprocessing.py +++ b/qsdsan/sanunits/_hydroprocessing.py @@ -50,18 +50,21 @@ class Hydrocracking(Reactor): HC catalyst lifetime, [hr]. hydrogen_P: float Hydrogen pressure, [Pa]. - hydrogen_rxned_to_heavy_oil: float - Reacted H2 to heavy oil mass ratio. + hydrogen_rxned_to_inf_oil: float + Reacted H2 to influent oil mass ratio. hydrogen_excess: float Actual hydrogen amount = hydrogen_rxned_to_biocrude*hydrogen_excess - hydrocarbon_ratio: float - Mass ratio of produced hydrocarbon to the sum of heavy oil and reacted H2. + oil_yield: float + Mass ratio of cracked oil to the sum of heavy oil and reacted H2, + gas yield is calculated as 1-oil_yield (about 100% conversion as in [1]). HCin_T: float HC influent temperature, [K]. HCrxn_T: float HC effluent (after reaction) temperature, [K]. - HC_composition: dict - HC effluent composition. + gas_composition: dict + Composition of the gas products, will be normalized to 100% sum. + oil_composition: dict + Composition of the cracked oil, will be normalized to 100% sum. References ---------- @@ -86,23 +89,22 @@ def __init__(self, ID='', ins=None, outs=(), thermo=None, WHSV=0.625, # wt./hr per wt. catalyst [1] catalyst_lifetime=5*7920, # 5 years [1] hydrogen_P=1039.7*6894.76, - hydrogen_rxned_to_heavy_oil=0.01125, + hydrogen_rxned_to_inf_oil=0.01125, hydrogen_excess=5.556, - hydrocarbon_ratio=1, # 100 wt% of heavy oil and reacted H2 - # nearly all input heavy oils and H2 will be converted to - # products [1] - # spreadsheet HC calculation + oil_yield=1-0.03880-0.00630, HCin_T=394+273.15, HCrxn_T=451+273.15, - HC_composition={'CO2':0.03880, 'CH4':0.00630, - 'CYCHEX':0.03714, 'HEXANE':0.01111, - 'HEPTANE':0.11474, 'OCTANE':0.08125, - 'C9H20':0.09086, 'C10H22':0.11756, - 'C11H24':0.16846, 'C12H26':0.13198, - 'C13H28':0.09302, 'C14H30':0.04643, - 'C15H32':0.03250, 'C16H34':0.01923, - 'C17H36':0.00431, 'C18H38':0.00099, - 'C19H40':0.00497, 'C20H42':0.00033}, + gas_composition={'CO2':0.03880, 'CH4':0.00630,}, + oil_composition={ + 'CYCHEX':0.03714, 'HEXANE':0.01111, + 'HEPTANE':0.11474, 'OCTANE':0.08125, + 'C9H20':0.09086, 'C10H22':0.11756, + 'C11H24':0.16846, 'C12H26':0.13198, + 'C13H28':0.09302, 'C14H30':0.04643, + 'C15H32':0.03250, 'C16H34':0.01923, + 'C17H36':0.00431, 'C18H38':0.00099, + 'C19H40':0.00497, 'C20H42':0.00033, + }, #combine C20H42 and PHYTANE as C20H42 # will not be a variable in uncertainty/sensitivity analysis P=None, tau=5, void_fraciton=0.4, # Towler @@ -116,12 +118,14 @@ def __init__(self, ID='', ins=None, outs=(), thermo=None, self.WHSV = WHSV self.catalyst_lifetime = catalyst_lifetime self.hydrogen_P = hydrogen_P - self.hydrogen_rxned_to_heavy_oil = hydrogen_rxned_to_heavy_oil + self.hydrogen_rxned_to_inf_oil = hydrogen_rxned_to_inf_oil self.hydrogen_excess = hydrogen_excess - self.hydrocarbon_ratio = hydrocarbon_ratio + self.oil_yield = oil_yield self.HCin_T = HCin_T + self._mixed_in = Stream(f'{ID}_mixed_in') self.HCrxn_T = HCrxn_T - self.HC_composition = HC_composition + self.gas_composition = gas_composition + self.oil_composition = oil_composition IC_in = Stream(f'{ID}_IC_in') IC_out = Stream(f'{ID}_IC_out') self.compressor = IsothermalCompressor(ID=f'.{ID}_IC', ins=IC_in, @@ -157,19 +161,20 @@ def _run(self): # catalysts amount is quite low compared to the main stream, therefore do not consider # heating/cooling of catalysts - hydrogen.imass['H2'] = heavy_oil.F_mass*self.hydrogen_rxned_to_heavy_oil*self.hydrogen_excess + hydrogen_rxned_to_inf_oil = self.hydrogen_rxned_to_inf_oil + hydrogen.imass['H2'] = heavy_oil.F_mass*hydrogen_rxned_to_inf_oil*self.hydrogen_excess hydrogen.phase = 'g' - hydrocarbon_mass = heavy_oil.F_mass*(1 +\ - self.hydrogen_rxned_to_heavy_oil)*\ - self.hydrocarbon_ratio - + hydrocarbon_mass = heavy_oil.F_mass*(1 + hydrogen_rxned_to_inf_oil) + # 100 wt% of heavy oil and reacted H2 + # nearly all input heavy oils and H2 will be converted to products [1] + # spreadsheet HC calculation hc_out.phase = 'g' - + for name, ratio in self.HC_composition.items(): hc_out.imass[name] = hydrocarbon_mass*ratio - hc_out.imass['H2'] = heavy_oil.F_mass*self.hydrogen_rxned_to_heavy_oil*(self.hydrogen_excess - 1) + hc_out.imass['H2'] = heavy_oil.F_mass*hydrogen_rxned_to_inf_oil*(self.hydrogen_excess - 1) hc_out.P = heavy_oil.P hc_out.T = self.HCrxn_T @@ -189,6 +194,36 @@ def _run(self): # make sure that carbon mass balance is within +/- 5%. Otherwise, an # exception will be raised. + def _normalize_composition(self, dct): + total = sum(dct.values()) + if total <=0: raise ValueError(f'Sum of total yields/composition should be positive, not {total}.') + return {k:v/total for k, v in dct.items()} + + @property + def gas_composition(self): + return self._gas_composition + @gas_composition.setter + def gas_composition(self, comp_dct): + self._gas_composition = self._normalize_composition(comp_dct) + + @property + def oil_composition(self): + return self._oil_composition + @oil_composition.setter + def oil_composition(self, comp_dct): + self._oil_composition = self._normalize_composition(comp_dct) + + @property + def HC_composition(self): + '''Composition of gas and oil products, normalized to 100%.''' + gas_composition = self.gas_composition + oil_composition = self.oil_composition + oil_yield = self.oil_yield + gas_yield = 1 - oil_yield + HC_composition = {k:v*gas_yield for k, v in gas_composition.items()} + HC_composition.update({k:v*oil_yield for k, v in oil_composition.items()}) + return self._normalize_composition(HC_composition) + @property def hydrocarbon_C(self): return sum(self.outs[0].imass[self.HC_composition]* @@ -208,6 +243,8 @@ def _design(self): hx_H2_ins0.copy_like(self.ins[1]) hx_H2_outs0.copy_like(hx_H2_ins0) hx_H2_ins0.phase = hx_H2_outs0.phase = 'g' + self._mixed_in.mix_from(self.ins) + if not self.HCin_T: self.HCin_T = self._mixed_in.T hx_H2_outs0.T = self.HCin_T hx_H2_ins0.P = hx_H2_outs0.P = IC_outs0.P hx_H2.simulate_as_auxiliary_exchanger(ins=hx_H2.ins, outs=hx_H2.outs) @@ -257,8 +294,8 @@ class Hydrotreating(Reactor): HT catalyst lifetime, [hr]. hydrogen_P: float Hydrogen pressure, [Pa]. - hydrogen_rxned_to_biocrude: float - Reacted H2 to biocrude mass ratio. + hydrogen_rxned_to_inf_oil: float + Reacted H2 to influent oil mass ratio. hydrogen_excess: float Actual hydrogen amount = hydrogen_rxned_to_biocrude*hydrogen_excess hydrocarbon_ratio: float @@ -303,7 +340,7 @@ def __init__(self, ID='', ins=None, outs=(), thermo=None, WHSV=0.625, # wt./hr per wt. catalyst [1] catalyst_lifetime=2*7920, # 2 years [1] hydrogen_P=1530*6894.76, - hydrogen_rxned_to_biocrude=0.046, + hydrogen_rxned_to_inf_oil=0.046, hydrogen_excess=3, hydrocarbon_ratio=0.875, # 87.5 wt% of biocrude and reacted H2 [1] # spreadsheet HT calculation @@ -354,7 +391,7 @@ def __init__(self, ID='', ins=None, outs=(), thermo=None, self.WHSV = WHSV self.catalyst_lifetime = catalyst_lifetime self.hydrogen_P = hydrogen_P - self.hydrogen_rxned_to_biocrude = hydrogen_rxned_to_biocrude + self.hydrogen_rxned_to_inf_oil = hydrogen_rxned_to_inf_oil self.hydrogen_excess = hydrogen_excess self.hydrocarbon_ratio = hydrocarbon_ratio self.HTin_T = HTin_T @@ -409,13 +446,13 @@ def _run(self): # heating/cooling of catalysts hydrogen_excess = self.hydrogen_excess - H2_rxned = biocrude.imass['Biocrude']*self.hydrogen_rxned_to_biocrude + H2_rxned = biocrude.imass['Biocrude']*self.hydrogen_rxned_to_inf_oil recovered_frac = (hydrogen_excess - 1)*self.PSA_efficiency*float(self.include_PSA) hydrogen.imass['H2'] = H2_rxned*(hydrogen_excess - recovered_frac) hydrogen.phase = 'g' hydrocarbon_mass = biocrude.imass['Biocrude']*\ - (1 + self.hydrogen_rxned_to_biocrude)*\ + (1 + self.hydrogen_rxned_to_inf_oil)*\ self.hydrocarbon_ratio ht_out.phase = 'g' From 0164999aab49f17153ade98c3ba78d87e20fd430 Mon Sep 17 00:00:00 2001 From: Yalin Date: Fri, 18 Oct 2024 11:12:01 -0400 Subject: [PATCH 12/32] more flexible catalyst ID in hydroprocessing units --- qsdsan/sanunits/_hydroprocessing.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/qsdsan/sanunits/_hydroprocessing.py b/qsdsan/sanunits/_hydroprocessing.py index ce89cef4..c075fbbb 100644 --- a/qsdsan/sanunits/_hydroprocessing.py +++ b/qsdsan/sanunits/_hydroprocessing.py @@ -48,6 +48,8 @@ class Hydrocracking(Reactor): Weight Hourly Space velocity, [kg feed/hr/kg catalyst]. catalyst_lifetime: float HC catalyst lifetime, [hr]. + catalyst_ID : str + ID of the catalyst. hydrogen_P: float Hydrogen pressure, [Pa]. hydrogen_rxned_to_inf_oil: float @@ -86,8 +88,10 @@ class Hydrocracking(Reactor): def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='Stream', + include_construction=False, WHSV=0.625, # wt./hr per wt. catalyst [1] catalyst_lifetime=5*7920, # 5 years [1] + catalyst_ID='HC_catalyst', hydrogen_P=1039.7*6894.76, hydrogen_rxned_to_inf_oil=0.01125, hydrogen_excess=5.556, @@ -114,9 +118,10 @@ def __init__(self, ID='', ins=None, outs=(), thermo=None, vessel_material='Stainless steel 316', vessel_type='Vertical'): - SanUnit.__init__(self, ID, ins, outs, thermo, init_with) + SanUnit.__init__(self, ID, ins, outs, thermo, init_with, include_construction=include_construction) self.WHSV = WHSV self.catalyst_lifetime = catalyst_lifetime + self.catalyst_ID = catalyst_ID self.hydrogen_P = hydrogen_P self.hydrogen_rxned_to_inf_oil = hydrogen_rxned_to_inf_oil self.hydrogen_excess = hydrogen_excess @@ -155,7 +160,7 @@ def _run(self): heavy_oil, hydrogen, catalyst_in = self.ins hc_out, catalyst_out = self.outs - catalyst_in.imass['HC_catalyst'] = heavy_oil.F_mass/self.WHSV/self.catalyst_lifetime + catalyst_in.imass[self.catalyst_ID] = heavy_oil.F_mass/self.WHSV/self.catalyst_lifetime catalyst_in.phase = 's' catalyst_out.copy_like(catalyst_in) # catalysts amount is quite low compared to the main stream, therefore do not consider @@ -292,6 +297,8 @@ class Hydrotreating(Reactor): Weight Hourly Space velocity, [kg feed/hr/kg catalyst]. catalyst_lifetime: float HT catalyst lifetime, [hr]. + catalyst_ID : str + ID of the catalyst. hydrogen_P: float Hydrogen pressure, [Pa]. hydrogen_rxned_to_inf_oil: float @@ -339,6 +346,7 @@ def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='Stream', WHSV=0.625, # wt./hr per wt. catalyst [1] catalyst_lifetime=2*7920, # 2 years [1] + catalyst_ID='HT_catalyst', hydrogen_P=1530*6894.76, hydrogen_rxned_to_inf_oil=0.046, hydrogen_excess=3, @@ -390,6 +398,7 @@ def __init__(self, ID='', ins=None, outs=(), thermo=None, SanUnit.__init__(self, ID, ins, outs, thermo, init_with) self.WHSV = WHSV self.catalyst_lifetime = catalyst_lifetime + self.catalyst_ID = catalyst_ID self.hydrogen_P = hydrogen_P self.hydrogen_rxned_to_inf_oil = hydrogen_rxned_to_inf_oil self.hydrogen_excess = hydrogen_excess @@ -439,7 +448,7 @@ def _run(self): HT_composition[chemical] /= (1-remove) HT_composition['PIPERDIN'] = 0 - catalyst_in.imass['HT_catalyst'] = biocrude.F_mass/self.WHSV/self.catalyst_lifetime + catalyst_in.imass[self.catalyst_ID] = biocrude.F_mass/self.WHSV/self.catalyst_lifetime catalyst_in.phase = 's' catalyst_out.copy_like(catalyst_in) # catalysts amount is quite low compared to the main stream, therefore do not consider From 5e2d8f0ec570628dc88d1d67b91f4fb791fe4c67 Mon Sep 17 00:00:00 2001 From: Yalin Date: Mon, 21 Oct 2024 07:54:04 -0400 Subject: [PATCH 13/32] improve setting in CHP utility stream prices --- qsdsan/sanunits/_combustion.py | 31 +++++++++++++++++++++++++++++-- 1 file changed, 29 insertions(+), 2 deletions(-) diff --git a/qsdsan/sanunits/_combustion.py b/qsdsan/sanunits/_combustion.py index 50602b8d..c5276a1a 100644 --- a/qsdsan/sanunits/_combustion.py +++ b/qsdsan/sanunits/_combustion.py @@ -21,6 +21,7 @@ # %% +import biosteam as bst from warnings import warn from flexsolve import IQ_interpolation from biosteam import HeatUtility, Facility @@ -302,7 +303,10 @@ def _cost(self): unit_CAPEX = self.unit_CAPEX unit_CAPEX /= 3600 # convert to $ per kJ/hr self.baseline_purchase_costs['CHP'] = unit_CAPEX * self.H_net_feeds - + # Update biosteam utility costs + uprices = bst.stream_utility_prices + uprices['Fuel'] = uprices['Natural gas'] = self.ins[1].price + uprices['Ash disposal'] = self.outs[1].price def _refresh_sys(self): sys = self._system @@ -311,10 +315,33 @@ def _refresh_sys(self): hu_dct = self._sys_heating_utilities = {} pu_dct = self._sys_power_utilities = {} for u in units: - hu_dct[u.ID] = tuple([i for i in u.heat_utilities if i.duty*i.flow>0]) + hu_dct[u.ID] = tuple([i for i in u.heat_utilities if ( + i.duty*i.flow>0 and i.agent.F_mass==i.agent.imass['H2O'])] + ) pu_dct[u.ID] = u.power_utility + @property + def fuel_price(self): + ''' + [Float] Price of fuel (natural gas), set to be the same as the price of ins[1] + and `bst.stream_utility_prices['Natural gas']`. + ''' + return self.ins[1].price + + natural_gas_price = fuel_price + + @property + def ash_disposal_price(self): + ''' + [Float] Price of ash disposal, set to be the same as the price of outs[1] + and `bst.stream_utility_prices['Ash disposal']`. + Negative means need to pay for ash disposal. + ''' + + """[Float] Price of ash disposal, same as `bst.stream_utility_prices['Ash disposal']`.""" + return self.outs[1].price + @property def CHP_type(self): ''' From 0613ba9d9f2876541a3c0495f7cf97a4e93d2011 Mon Sep 17 00:00:00 2001 From: Yalin Date: Mon, 21 Oct 2024 08:33:29 -0400 Subject: [PATCH 14/32] fix minor doc issue --- docs/source/FAQ.rst | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/docs/source/FAQ.rst b/docs/source/FAQ.rst index 131a847d..54e9ab30 100644 --- a/docs/source/FAQ.rst +++ b/docs/source/FAQ.rst @@ -77,6 +77,28 @@ There are multiple possible reasons: Then when you open the Jupyter Notebook, select the ```` kernel when you create a new notebook you can find more details in this post about `enabling multiple kernels in Jupyter Notebook `_. +``Underlying object has vanished`` +********************************** +This error is related to ``numba`` caching, we haven't figured out the exact mechanism, but clearing cache will help resolve it. One/both of the following approaches should work: + +1. Clear cache. Remove all ``.pyc``, ``.nbc``, and ``.nbi`` files, you can do this in your CLI using (replace with the directory to your ``thermosteam``, ``biosteam``, ``qsdsan``, and ``exposan`` directory): + + .. code:: + + get-childitem -recurse -include *.pyc | remove-item + get-childitem -recurse -include *.nbc | remove-item + get-childitem -recurse -include *.nbi | remove-item + +2. Uninstalling and reinstalling a different version of ``numba``. Suppose you now have 0.58.1 and the newest version is 0.60.0, you can do: + + .. code:: + + pip uninstall numba + pip install --no-cache-dir numba==0.60.0 + +The ``--no-cache-dir`` option is to do a fresh installation rather than using previously downloaded packages. Note that you need to exit out your editor/any other programs that are currently using numba. Otherwise the uninstallation is incomplete, you might be prompted to do a manual removal, or this won't work. + + ``UnicodeDecodeError`` ********************** When using non-English operating systems, you may run into errors similar to (cp949 is the case of Korean Windows): From 0be6a66eef2ebc0d78633dca4a9e67d1c6ee3d1f Mon Sep 17 00:00:00 2001 From: Yalin Date: Mon, 21 Oct 2024 08:34:10 -0400 Subject: [PATCH 15/32] fix minor doc issue --- qsdsan/sanunits/_reactor.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/qsdsan/sanunits/_reactor.py b/qsdsan/sanunits/_reactor.py index 89d5cec1..029daac4 100644 --- a/qsdsan/sanunits/_reactor.py +++ b/qsdsan/sanunits/_reactor.py @@ -57,7 +57,7 @@ class Reactor(SanUnit, PressureVessel, isabstract=True): Power usage of agitator (converted from 0.5 hp/1000 gal as in [1]). If mixing_intensity is provided, this will be calculated based on - the mixing_intensity and viscosity of the influent mixture as in [2]_ + the mixing_intensity and viscosity of the influent mixture as in [2] wall_thickness_factor=1: float A safety factor to scale up the calculated minimum wall thickness. vessel_material : str, optional @@ -67,10 +67,11 @@ class Reactor(SanUnit, PressureVessel, isabstract=True): References ---------- - .. [1] Seider, W. D.; Lewin, D. R.; Seader, J. D.; Widagdo, S.; Gani, R.; + [1] Seider, W. D.; Lewin, D. R.; Seader, J. D.; Widagdo, S.; Gani, R.; Ng, M. K. Cost Accounting and Capital Cost Estimation. In Product and Process Design Principles; Wiley, 2017; pp 470. - .. [2] Shoener et al. Energy Positive Domestic Wastewater Treatment: + + [2] Shoener et al. Energy Positive Domestic Wastewater Treatment: The Roles of Anaerobic and Phototrophic Technologies. Environ. Sci.: Processes Impacts 2014, 16 (6), 1204–1222. https://doi.org/10.1039/C3EM00711A. From 8914bfc8965fddd0229e53af8ea4083f452296dd Mon Sep 17 00:00:00 2001 From: Yalin Date: Tue, 22 Oct 2024 10:37:02 -0400 Subject: [PATCH 16/32] better steam vs. natural gas utility accounting in `CHP` --- qsdsan/sanunits/_combustion.py | 51 ++++++++++++++++++++++++---------- qsdsan/utils/utilities.py | 10 +++++++ 2 files changed, 46 insertions(+), 15 deletions(-) diff --git a/qsdsan/sanunits/_combustion.py b/qsdsan/sanunits/_combustion.py index c5276a1a..f9b6dc23 100644 --- a/qsdsan/sanunits/_combustion.py +++ b/qsdsan/sanunits/_combustion.py @@ -185,6 +185,7 @@ def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='WasteStream self.system = None self.supplement_power_utility = supplement_power_utility self._sys_heating_utilities = () + self._sys_steam_utilities = () self._sys_power_utilities = () def _init_lca(self): @@ -239,11 +240,11 @@ def react(natural_gas_flow=0): # Calculate all energy needs in kJ/hr as in H_net_feeds kwds = dict(system=self.system, operating_hours=self.system.operating_hours, exclude_units=(self,)) pu = self.power_utility - H_heating_needs = sum_system_utility(**kwds, utility='heating', result_unit='kJ/hr')/self.combustion_eff + H_steam_needs = sum_system_utility(**kwds, utility='steam', result_unit='kJ/hr')/self.combustion_eff H_power_needs = sum_system_utility(**kwds, utility='power', result_unit='kJ/hr')/self.combined_eff # Calculate the amount of energy needs to be provided - H_supp = H_heating_needs+H_power_needs if self.supplement_power_utility else H_heating_needs + H_supp = H_steam_needs+H_power_needs if self.supplement_power_utility else H_steam_needs # Objective function to calculate the heat deficit at a given natural gas flow rate def H_deficit_at_natural_gas_flow(flow): @@ -264,17 +265,19 @@ def H_deficit_at_natural_gas_flow(flow): H_net_feeds = react(0) # Update heating utilities - self.heat_utilities = HeatUtility.sum_by_agent(sum(self.sys_heating_utilities.values(), ())) + self.steam_utilities = HeatUtility.sum_by_agent(sum(self.sys_steam_utilities.values(), [])) + ngu = self.natural_gas_utilities = HeatUtility.sum_by_agent(sum(self.sys_natural_gas_utilities.values(), [])) + self.heat_utilities = HeatUtility.sum_by_agent(sum(self.sys_heating_utilities.values(), [])) + natural_gas.imol['CH4'] += sum(i.flow for i in ngu) # natural gas is added on separately for hu in self.heat_utilities: hu.reverse() - # Power production if there is sufficient energy - if H_net_feeds <= H_heating_needs: + if H_net_feeds <= H_steam_needs: pu.production = 0 else: - pu.production = (H_net_feeds-H_heating_needs)/3600*self.combined_eff + pu.production = (H_net_feeds-H_steam_needs)/3600*self.combined_eff - self.H_heating_needs = H_heating_needs + self.H_steam_needs = H_steam_needs self.H_power_needs = H_power_needs self.H_net_feeds = H_net_feeds @@ -310,16 +313,24 @@ def _cost(self): def _refresh_sys(self): sys = self._system + ng_dct = self._sys_natural_gas_utilities = {} + steam_dct = self._sys_steam_utilities = {} + pu_dct = self._sys_power_utilities = {} if sys: units = [u for u in sys.units if u is not self] - hu_dct = self._sys_heating_utilities = {} - pu_dct = self._sys_power_utilities = {} for u in units: - hu_dct[u.ID] = tuple([i for i in u.heat_utilities if ( - i.duty*i.flow>0 and i.agent.F_mass==i.agent.imass['H2O'])] - ) - pu_dct[u.ID] = u.power_utility - + pu = u.power_utility + if pu: pu_dct[u.ID] = pu + steam_utilities = [] + for hu in u.heat_utilities: + if hu.flow*hu.duty <= 0: continue # cooling utilities + if hu.ID=='natural_gas': ng_dct[u.ID] = [hu] + elif 'steam' in hu.ID: steam_utilities.append(hu) + else: raise ValueError(f'The heating utility {hu.ID} is not recognized by the CHP.') + if steam_utilities: steam_dct[u.ID] = steam_utilities + sys_hus = {k:ng_dct.get(k,[])+steam_dct.get(k, []) + for k in list(ng_dct.keys())+list(steam_dct.keys())} + self._sys_heating_utilities = sys_hus @property def fuel_price(self): @@ -427,9 +438,19 @@ def system(self, i): @property def sys_heating_utilities(self): - '''[dict] Heating utilities of the given system (excluding this CHP unit).''' + '''[dict] Heating utilities (steams and natural gas) of the given system (excluding this CHP unit).''' return self._sys_heating_utilities + @property + def sys_natural_gas_utilities(self): + '''[dict] Steam utilities of the given system (excluding this CHP unit).''' + return self._sys_natural_gas_utilities + + @property + def sys_steam_utilities(self): + '''[dict] Steam utilities of the given system (excluding this CHP unit).''' + return self._sys_steam_utilities + @property def sys_power_utilities(self): '''[dict] Power utilities of the given system (excluding this CHP unit).''' diff --git a/qsdsan/utils/utilities.py b/qsdsan/utils/utilities.py index b8c61791..da92a43c 100644 --- a/qsdsan/utils/utilities.py +++ b/qsdsan/utils/utilities.py @@ -48,6 +48,8 @@ def sum_system_utility(system, operating_hours=None, exclude_units=(), >>> from qsdsan.utils import create_example_system, sum_system_utility >>> sys = create_example_system() >>> sys.simulate() + >>> sum_system_utility(sys, utility='steam', result_unit='kJ/yr') # doctest: +ELLIPSIS + 463479... >>> sum_system_utility(sys, utility='heating', result_unit='kJ/yr') # doctest: +ELLIPSIS 463479... >>> sum_system_utility(sys, utility='cooling', result_unit='GJ/yr') # doctest: +NUMBER @@ -70,6 +72,14 @@ def sum_system_utility(system, operating_hours=None, exclude_units=(), attr = 'consumption' if not calculate_net_utility else 'rate' tot = sum([get(i.power_utility, attr) for i in units if i.power_utility])*hrs return auom('kWh/yr').convert(tot, unit) + elif utility == 'steam': + unit = 'GJ/yr' if not result_unit else result_unit + hu = sum([i.heat_utilities for i in units], []) + if not calculate_net_utility: + tot = sum([i.duty for i in hu if 'steam' in i.ID])/1e6*hrs + else: + tot = sum([i.duty for i in hu if ('steam' in i.ID and i.duty>0)])/1e6*hrs + return auom('GJ/yr').convert(tot, unit) elif utility == 'heating': unit = 'GJ/yr' if not result_unit else result_unit hu = sum([i.heat_utilities for i in units], []) From bce6a2b925bfe7b85fcf08db2ecfe411fcca0f84 Mon Sep 17 00:00:00 2001 From: Yalin Date: Tue, 22 Oct 2024 23:17:42 -0400 Subject: [PATCH 17/32] fix minor bug in `Reactor` --- qsdsan/sanunits/_reactor.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/qsdsan/sanunits/_reactor.py b/qsdsan/sanunits/_reactor.py index 029daac4..79a6c67b 100644 --- a/qsdsan/sanunits/_reactor.py +++ b/qsdsan/sanunits/_reactor.py @@ -218,7 +218,7 @@ def vessel_material(self, i): exist_material = getattr(self, '_vessel_material', None) PressureVessel.vessel_material.fset(self, i) if i and exist_material == i: return # type doesn't change, no need to reload construction items - self._init_lca() + if self.include_construction: self._init_lca() @property def kW_per_m3(self): @@ -230,7 +230,6 @@ def kW_per_m3(self): mixture.mix_from(self.ins) kW_per_m3 = mixture.mu*(G**2)/1e3 return kW_per_m3 - @kW_per_m3.setter def kW_per_m3(self, i): if self.mixing_intensity and i is not None: From 066f42f527715a2671243d59cfc1697a9d6b2d29 Mon Sep 17 00:00:00 2001 From: Yalin Date: Thu, 24 Oct 2024 08:46:14 -0400 Subject: [PATCH 18/32] fix bug for biosteam units without results --- qsdsan/_sanunit.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/qsdsan/_sanunit.py b/qsdsan/_sanunit.py index 23a5adfb..51cac1a7 100644 --- a/qsdsan/_sanunit.py +++ b/qsdsan/_sanunit.py @@ -631,6 +631,7 @@ def results(self, with_units=True, include_utilities=True, include_total_cost=True, include_installed_cost=False, include_zeros=True, external_utilities=(), key_hook=None): + if super().results is None: return super().results results = super().results(with_units, include_utilities, include_total_cost, include_installed_cost, include_zeros, external_utilities, key_hook) @@ -647,7 +648,7 @@ def results(self, with_units=True, include_utilities=True, results.insert(0, 'Units', '') results.loc[(k, ''), :] = ('USD/hr', v) results.columns.name = type(self).__name__ - if with_units: + if with_units and results is not None: results.replace({'USD': f'{currency}', 'USD/hr': f'{currency}/hr'}, inplace=True) return results From 6bb5b855068946f094b75f5e47537ff7c82c1264 Mon Sep 17 00:00:00 2001 From: Yalin Date: Thu, 31 Oct 2024 06:20:40 -0700 Subject: [PATCH 19/32] add TEA indices --- qsdsan/__init__.py | 2 +- qsdsan/utils/__init__.py | 3 + qsdsan/utils/indices.py | 198 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 202 insertions(+), 1 deletion(-) create mode 100644 qsdsan/utils/indices.py diff --git a/qsdsan/__init__.py b/qsdsan/__init__.py index 276f3f79..c667a18c 100644 --- a/qsdsan/__init__.py +++ b/qsdsan/__init__.py @@ -45,7 +45,6 @@ Flowsheet = _bst.Flowsheet main_flowsheet = _bst.main_flowsheet default_utilities = _bst.default_utilities -CEPCI_by_year = _bst.units.design_tools.CEPCI_by_year # Global variables currency = 'USD' @@ -54,6 +53,7 @@ from . import utils +CEPCI_by_year = utils.indices.tea_indices['CEPCI'] from ._component import * from ._components import * from ._sanstream import * diff --git a/qsdsan/utils/__init__.py b/qsdsan/utils/__init__.py index 5f5317b3..af8be569 100644 --- a/qsdsan/utils/__init__.py +++ b/qsdsan/utils/__init__.py @@ -47,6 +47,7 @@ formatting, loading, dynamics, + indices, misc, model_eval, parsing, @@ -62,6 +63,7 @@ from .formatting import * from .loading import * from .dynamics import * +from .indices import * from .misc import * from .model_eval import * from .parsing import * @@ -79,6 +81,7 @@ *formatting.__all__, *loading.__all__, *dynamics.__all__, + *indices.__all__, *model_eval.__all__, *misc.__all__, *parsing.__all__, diff --git a/qsdsan/utils/indices.py b/qsdsan/utils/indices.py new file mode 100644 index 00000000..f5b7afad --- /dev/null +++ b/qsdsan/utils/indices.py @@ -0,0 +1,198 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +''' +QSDsan: Quantitative Sustainable Design for sanitation and resource recovery systems + +This module is developed by: + + Yalin Li + +This module is under the University of Illinois/NCSA Open Source License. +Please refer to https://github.com/QSD-Group/QSDsan/blob/main/LICENSE.txt +for license details. +''' + +__all__ = ( + 'tea_indices', + ) + + +# %% + +CEPCI = { + 1990: 357.6, + 1991: 361.3, + 1992: 358.2, + 1993: 359.2, + 1994: 368.1, + 1995: 381.1, + 1996: 381.7, + 1997: 386.5, + 1998: 389.5, + 1999: 390.6, + 2000: 394.1, + 2001: 394.3, + 2002: 395.6, + 2003: 402.0, + 2004: 444.2, + 2005: 468.2, + 2006: 499.6, + 2007: 525.4, + 2008: 575.4, + 2009: 521.9, + 2010: 550.8, + 2011: 585.7, + 2012: 584.6, + 2013: 567.3, + 2014: 576.1, + 2015: 556.8, + 2016: 541.7, + 2017: 567.5, + 2018: 603.1, + 2019: 607.5, + 2020: 596.2, + 2021: 708.8, + 2022: 816.0, + 2023: 798.0, + } + +ChemPPI = { + 1984: 100.0, + 1985: 100.8, + 1986: 100.6, + 1987: 103.6, + 1988: 113.0, + 1989: 119.6, + 1990: 121.0, + 1991: 124.4, + 1992: 125.8, + 1993: 127.2, + 1994: 130.0, + 1995: 143.4, + 1996: 145.8, + 1997: 147.1, + 1998: 148.7, + 1999: 149.7, + 2000: 156.7, + 2001: 158.4, + 2002: 157.3, + 2003: 164.6, + 2004: 172.8, + 2005: 187.2, + 2006: 196.8, + 2007: 203.3, + 2008: 228.2, + 2009: 224.7, + 2010: 233.7, + 2011: 252.1, + 2012: 260.3, + 2013: 263.9, + 2014: 269.2, + 2015: 264.8, + 2016: 267.1, + 2017: 277.6, + 2018: 290.2, + 2019: 292.1, + 2020: 289.0, + 2021: 708.8, + 2022: 816.0, + 'Seider': 567.0, + } + +# U.S. Energy Information Administration (EIA) Annual Energy Outlook (AEO) +GDP = { + 2003: 0.808, + 2005: 0.867, + 2007: 0.913, + 2008: 0.941, + 2009: 0.951, + 2010: 0.962, + 2011: 0.983, + 2012: 1.000, + 2013: 1.014, + 2014: 1.033, + 2015: 1.046, + 2016: 1.059, + 2017: 1.078, + 2018: 1.100, + 2019: 1.123, + 2020: 1.133, + 2021: 1.181, + 2022: 1.269, + 2023: 1.322, + 2024: 1.354, + } + +# https://data.bls.gov/cgi-bin/srgate; CEU3232500008; Chemicals +labor = { + 1990: 12.85, + 1991: 13.30, + 1992: 13.70, + 1993: 13.97, + 1994: 14.33, + 1995: 14.86, + 1996: 15.37, + 1997: 15.78, + 1998: 16.23, + 1999: 16.40, + 2000: 17.09, + 2001: 17.57, + 2002: 17.97, + 2003: 18.50, + 2004: 19.17, + 2005: 19.67, + 2006: 19.60, + 2007: 19.55, + 2008: 19.50, + 2009: 20.30, + 2010: 21.07, + 2011: 21.45, + 2012: 21.45, + 2013: 21.40, + 2014: 21.49, + 2015: 21.76, + 2016: 22.72, + 2017: 24.28, + 2018: 25.46, + 2019: 25.46, + 2020: 26.04, + 2021: 26.69, + 2022: 27.35, + 2023: 29.77, + } + +# Federal Reserve Economic Data, Personal Consumption Expenditures: Chain-type Price Index, Index 2017=1.00, Annual, Seasonally Adjusted +PCEPI = { + 2000: 73.822, + 2001: 75.302, + 2002: 76.291, + 2003: 77.894 , + 2004: 79.827, + 2005: 82.127, + 2006: 84.440, + 2007: 86.607, + 2008: 89.170, + 2009: 88.921, + 2010: 90.514, + 2011: 92.804, + 2012: 94.534, + 2013: 95.781, + 2014: 97.121, + 2015: 97.299, + 2016: 98.284, + 2017: 100.000, + 2018: 102.047, + 2019: 103.513, + 2020: 104.635, + 2021: 109.001, + 2022: 116.043, + 2023: 120.380, + 2024: 121.966, + } + +tea_indices = { + 'CEPCI': CEPCI, + 'ChemPPI': ChemPPI, + 'labor': labor, + 'PCEPI': PCEPI, + } \ No newline at end of file From 95f17c3801cb96de2e8c28a3235f317e356b0202 Mon Sep 17 00:00:00 2001 From: Yalin Date: Sun, 3 Nov 2024 07:18:32 -0500 Subject: [PATCH 20/32] add TEA-related indices --- qsdsan/utils/indices.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/qsdsan/utils/indices.py b/qsdsan/utils/indices.py index f5b7afad..4d013460 100644 --- a/qsdsan/utils/indices.py +++ b/qsdsan/utils/indices.py @@ -7,6 +7,8 @@ Yalin Li + Ali Ahmad + This module is under the University of Illinois/NCSA Open Source License. Please refer to https://github.com/QSD-Group/QSDsan/blob/main/LICENSE.txt for license details. From b494632cb4524e6e9ac84d14e572f61db57cfd18 Mon Sep 17 00:00:00 2001 From: Yalin Date: Tue, 5 Nov 2024 08:57:02 -0500 Subject: [PATCH 21/32] better solution for cleaning numba cache --- docs/source/FAQ.rst | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/docs/source/FAQ.rst b/docs/source/FAQ.rst index 54e9ab30..f92f22ec 100644 --- a/docs/source/FAQ.rst +++ b/docs/source/FAQ.rst @@ -77,7 +77,7 @@ There are multiple possible reasons: Then when you open the Jupyter Notebook, select the ```` kernel when you create a new notebook you can find more details in this post about `enabling multiple kernels in Jupyter Notebook `_. -``Underlying object has vanished`` +``underlying object has vanished`` ********************************** This error is related to ``numba`` caching, we haven't figured out the exact mechanism, but clearing cache will help resolve it. One/both of the following approaches should work: @@ -85,9 +85,7 @@ This error is related to ``numba`` caching, we haven't figured out the exact mec .. code:: - get-childitem -recurse -include *.pyc | remove-item - get-childitem -recurse -include *.nbc | remove-item - get-childitem -recurse -include *.nbi | remove-item + get-childitem . -recurse -include *.pyc, *.nbc, *.nbi | remove-item 2. Uninstalling and reinstalling a different version of ``numba``. Suppose you now have 0.58.1 and the newest version is 0.60.0, you can do: From cac89b3623a3ee227374754f253854a6c9dd9b37 Mon Sep 17 00:00:00 2001 From: Yalin Date: Tue, 12 Nov 2024 15:04:26 -0500 Subject: [PATCH 22/32] fix bug in `SanStream.impact_item` --- qsdsan/_impact_item.py | 2 +- qsdsan/_sanstream.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/qsdsan/_impact_item.py b/qsdsan/_impact_item.py index 395bfb1c..b7c83d84 100644 --- a/qsdsan/_impact_item.py +++ b/qsdsan/_impact_item.py @@ -822,7 +822,7 @@ def linked_stream(self, new_s): f'is replaced with {self.ID}.') else: warn(f'The original `StreamImpactItem` linked to stream {new_s.ID} ' - f'is replaced with upon the creation of a new stream.') + f'is replaced upon the creation of a new stream.') new_s._stream_impact_item = self self._linked_stream = new_s diff --git a/qsdsan/_sanstream.py b/qsdsan/_sanstream.py index d51e53b8..c2beccb8 100644 --- a/qsdsan/_sanstream.py +++ b/qsdsan/_sanstream.py @@ -170,12 +170,13 @@ def copy_flow(self, other, IDs=..., *, remove=False, exclude=False): -------- :func:`copy` for the differences between ``copy``, ``copy_like``, and ``copy_flow``. ''' + stream_impact_item = self.stream_impact_item Stream.copy_flow(self, other=other, IDs=IDs, remove=remove, exclude=exclude) if not isinstance(other, SanStream): return - self._stream_impact_item = None + self._stream_impact_item = stream_impact_item def flow_proxy(self, ID=None): From d9b4f4edd87f8dd5a45946b1c053f8b696e773c3 Mon Sep 17 00:00:00 2001 From: Yalin Date: Tue, 12 Nov 2024 15:21:25 -0500 Subject: [PATCH 23/32] minor fix in LCA table formatting --- qsdsan/_lca.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qsdsan/_lca.py b/qsdsan/_lca.py index 1c69a169..9faecd60 100644 --- a/qsdsan/_lca.py +++ b/qsdsan/_lca.py @@ -755,7 +755,7 @@ def get_impact_table(self, category, annual=False): item_dct[key] = [] for other_ID in self.other_items.keys(): other = self.other_items[other_ID]['item'] - item_dct['Other'].append(f'{other_ID} [{other.functional_unit}]') + item_dct['Other'].append(f'{other_ID}') quantity = self.other_items[other_ID]['quantity'] item_dct['Quantity'].append(quantity) for ind in self.indicators: From ceb143982134a11388b7479466f9f8dc2834edb7 Mon Sep 17 00:00:00 2001 From: Yalin Date: Fri, 15 Nov 2024 10:55:57 -0500 Subject: [PATCH 24/32] fix `LCA` annual results --- qsdsan/_lca.py | 66 ++++++++++++++++++++++++++++++-------------------- 1 file changed, 40 insertions(+), 26 deletions(-) diff --git a/qsdsan/_lca.py b/qsdsan/_lca.py index 9faecd60..2e2710d5 100644 --- a/qsdsan/_lca.py +++ b/qsdsan/_lca.py @@ -175,6 +175,9 @@ class LCA: >>> # Retrieve impacts associated with a specific indicator >>> lca.get_total_impacts()[GWP.ID] # doctest: +ELLIPSIS 349737809... + >>> # Annual results + >>> lca.get_total_impacts(annual=True)[GWP.ID] # doctest: +ELLIPSIS + 34973780... >>> # Or breakdowns of the different category >>> lca.get_impact_table('Construction') # doctest: +SKIP >>> # Below is for testing purpose, you do not need it @@ -637,12 +640,13 @@ def get_unit_impacts( return tot - def _append_cat_sum(self, cat_table, cat, tot): + def _append_cat_sum(self, cat_table, cat, tot, annual=False): num = len(cat_table) cat_table.loc[num] = '' # initiate a blank spot for value to be added later - + suffix = '/yr' if annual else '' + for i in self.indicators: - cat_table[f'{i.ID} [{i.unit}]'][num] = tot[i.ID] + cat_table[f'{i.ID} [{i.unit}{suffix}]'][num] = tot[i.ID] cat_table[f'Category {i.ID} Ratio'][num] = 1 if cat in ('construction', 'transportation'): @@ -662,17 +666,21 @@ def get_impact_table(self, category, annual=False): Parameters ---------- category : str - Can be 'construction', 'transportation', 'stream', or 'other'. + Can be 'Construction', 'Transportation', 'Stream', or 'Other'. annual : bool If True, will return the annual impacts considering `uptime_ratio` instead of across the system lifetime. ''' time = self.lifetime_hr + sys_yr = self.lifetime cat = category.lower() tot_f = getattr(self, f'get_{cat}_impacts') - kwargs = {'annual': annual} if cat != 'other' else {} + # kwargs = {'annual': annual} if cat != 'other' else {} + kwargs = {'annual': annual} tot = tot_f(**kwargs) - + suffix = '/yr' if annual else'' + _append_cat_sum = self._append_cat_sum + if cat in ('construction', 'transportation'): units = sorted(getattr(self, f'_{cat}_units'), key=(lambda su: su.ID)) @@ -684,31 +692,35 @@ def get_impact_table(self, category, annual=False): # Note that item_dct = dict.fromkeys([item.ID for item in items], []) won't work item_dct = dict.fromkeys([item.ID for item in items]) for item_ID in item_dct.keys(): - item_dct[item_ID] = dict(SanUnit=[], Quantity=[]) + item_dct[item_ID] = {'SanUnit': [], f'Quantity{suffix}': []} for su in units: if not isinstance(su, SanUnit): continue for i in getattr(su, cat): item_dct[i.item.ID]['SanUnit'].append(su.ID) if cat == 'transportation': - item_dct[i.item.ID]['Quantity'].append(i.quantity*time/i.interval) + quantity = i.quantity*time/i.interval + quantity = quantity/sys_yr if annual else quantity + item_dct[i.item.ID][f'Quantity{suffix}'].append(quantity) else: # construction lifetime = i.lifetime or su.lifetime or self.lifetime if isinstance(lifetime, dict): # in the case the the equipment is not in the unit lifetime dict lifetime = lifetime.get(i.item.ID) or self.lifetime - constr_ratio = self.lifetime/lifetime if self.annualize_construction else ceil(self.lifetime/lifetime) - item_dct[i.item.ID]['Quantity'].append(i.quantity*constr_ratio) + constr_ratio = sys_yr/lifetime if self.annualize_construction else ceil(sys_yr/lifetime) + quantity = i.quantity * constr_ratio + quantity = quantity/sys_yr if annual else quantity + item_dct[i.item.ID][f'Quantity{suffix}'].append(quantity) dfs = [] for item in items: dct = item_dct[item.ID] dct['SanUnit'].append('Total') - dct['Quantity'] = np.append(dct['Quantity'], sum(dct['Quantity'])) - if dct['Quantity'].sum() == 0.: dct['Item Ratio'] = 0 - else: dct['Item Ratio'] = dct['Quantity']/dct['Quantity'].sum()*2 + dct[f'Quantity{suffix}'] = np.append(dct[f'Quantity{suffix}'], sum(dct[f'Quantity{suffix}'])) + if dct[f'Quantity{suffix}'].sum() == 0.: dct['Item Ratio'] = 0 + else: dct['Item Ratio'] = dct[f'Quantity{suffix}']/dct[f'Quantity{suffix}'].sum()*2 for i in self.indicators: if i.ID in item.CFs: - dct[f'{i.ID} [{i.unit}]'] = impact = dct['Quantity']*item.CFs[i.ID] + dct[f'{i.ID} [{i.unit}{suffix}]'] = impact = dct[f'Quantity{suffix}']*item.CFs[i.ID] dct[f'Category {i.ID} Ratio'] = impact/tot[i.ID] else: dct[f'{i.ID} [{i.unit}]'] = dct[f'Category {i.ID} Ratio'] = 0 @@ -721,13 +733,13 @@ def get_impact_table(self, category, annual=False): dfs.append(df) table = pd.concat(dfs) - return self._append_cat_sum(table, cat, tot) + return _append_cat_sum(table, cat, tot, annual=annual) - ind_head = sum(([f'{i.ID} [{i.unit}]', + ind_head = sum(([f'{i.ID} [{i.unit}{suffix}]', f'Category {i.ID} Ratio'] for i in self.indicators), []) if cat in ('stream', 'streams'): - headings = ['Stream', 'Mass [kg]', *ind_head] + headings = ['Stream', f'Mass [kg]{suffix}', *ind_head] item_dct = dict.fromkeys(headings) for key in item_dct.keys(): item_dct[key] = [] @@ -735,21 +747,22 @@ def get_impact_table(self, category, annual=False): ws = ws_item.linked_stream item_dct['Stream'].append(ws.ID) mass = ws_item.flow_getter(ws) * time - item_dct['Mass [kg]'].append(mass) + mass = mass/sys_yr if annual else mass + item_dct[f'Mass [kg]{suffix}'].append(mass) for ind in self.indicators: if ind.ID in ws_item.CFs.keys(): impact = ws_item.CFs[ind.ID]*mass - item_dct[f'{ind.ID} [{ind.unit}]'].append(impact) + item_dct[f'{ind.ID} [{ind.unit}{suffix}]'].append(impact) item_dct[f'Category {ind.ID} Ratio'].append(impact/tot[ind.ID]) else: - item_dct[f'{ind.ID} [{ind.unit}]'].append(0) + item_dct[f'{ind.ID} [{ind.unit}{suffix}]'].append(0) item_dct[f'Category {ind.ID} Ratio'].append(0) table = pd.DataFrame.from_dict(item_dct) table.set_index(['Stream'], inplace=True) - return self._append_cat_sum(table, cat, tot) + return _append_cat_sum(table, cat, tot, annual=annual) elif cat == 'other': - headings = ['Other', 'Quantity', *ind_head] + headings = ['Other', f'Quantity{suffix}', *ind_head] item_dct = dict.fromkeys(headings) for key in item_dct.keys(): item_dct[key] = [] @@ -757,19 +770,20 @@ def get_impact_table(self, category, annual=False): other = self.other_items[other_ID]['item'] item_dct['Other'].append(f'{other_ID}') quantity = self.other_items[other_ID]['quantity'] - item_dct['Quantity'].append(quantity) + quantity = quantity/sys_yr if annual else quantity + item_dct[f'Quantity{suffix}'].append(quantity) for ind in self.indicators: if ind.ID in other.CFs.keys(): impact = other.CFs[ind.ID]*quantity - item_dct[f'{ind.ID} [{ind.unit}]'].append(impact) + item_dct[f'{ind.ID} [{ind.unit}{suffix}]'].append(impact) item_dct[f'Category {ind.ID} Ratio'].append(impact/tot[ind.ID]) else: - item_dct[f'{ind.ID} [{ind.unit}]'].append(0) + item_dct[f'{ind.ID} [{ind.unit}{suffix}]'].append(0) item_dct[f'Category {ind.ID} Ratio'].append(0) table = pd.DataFrame.from_dict(item_dct) table.set_index(['Other'], inplace=True) - return self._append_cat_sum(table, cat, tot) + return _append_cat_sum(table, cat, tot, annual=annual) raise ValueError( 'category can only be "Construction", "Transportation", "Stream", or "Other", ' \ From 19486fa18259392f1b23de9bb6d6b8d2ca906442 Mon Sep 17 00:00:00 2001 From: Yalin Date: Tue, 26 Nov 2024 15:21:20 -0500 Subject: [PATCH 25/32] more info on chemical indices --- qsdsan/utils/indices.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/qsdsan/utils/indices.py b/qsdsan/utils/indices.py index 4d013460..7c575a2e 100644 --- a/qsdsan/utils/indices.py +++ b/qsdsan/utils/indices.py @@ -125,7 +125,7 @@ 2024: 1.354, } -# https://data.bls.gov/cgi-bin/srgate; CEU3232500008; Chemicals +# https://data.bls.gov/cgi-bin/srgate; CEU3232500008; Chemical manufacturing labor = { 1990: 12.85, 1991: 13.30, @@ -168,7 +168,7 @@ 2000: 73.822, 2001: 75.302, 2002: 76.291, - 2003: 77.894 , + 2003: 77.894, 2004: 79.827, 2005: 82.127, 2006: 84.440, From 692c94c12887232ec6200d847327fc564e718c26 Mon Sep 17 00:00:00 2001 From: Yalin Li Date: Sat, 7 Dec 2024 17:38:22 -0500 Subject: [PATCH 26/32] fix invalid escape by adding "r" before docs, more details in https://github.com/HandBrake/HandBrake/issues/5454 --- qsdsan/stats.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/qsdsan/stats.py b/qsdsan/stats.py index 741f4b82..7c032891 100644 --- a/qsdsan/stats.py +++ b/qsdsan/stats.py @@ -405,7 +405,7 @@ def morris_till_convergence(model, inputs, metrics=None, nan_policy='propagate', conf_level=0.95, print_to_console=False, file='', **kwargs): - ''' + r''' Run Morris analysis from N=2 to N=N_max until the results converge (i.e., mu_star_conf/mu_star_max < threshold for all parameters, where as mu_star_max is the maximum :math:`{\mu^*}` value for a certain metric, @@ -1088,7 +1088,7 @@ def plot_morris_results(morris_dct, metric, kind='scatter', ax=None, x_axis='mu_star', plot_lines=True, k1=0.1, k2=0.5, k3=1, label_kind='number', color='k', file='', close_fig=True, **kwargs): - ''' + r''' Visualize the results from Morris One-at-A-Time analysis as either scatter or bar plots. In scatter plots, the x values are :math:`{\mu^*}` and the y values are :math:`{\sigma}`. @@ -1213,7 +1213,7 @@ def plot_morris_results(morris_dct, metric, kind='scatter', ax=None, def plot_morris_convergence(result_dct, metric, parameters=(), plot_rank=False, kind='line', ax=None, show_error=True, palette='pastel', file='', close_fig=True): - ''' + r''' Plot the evolution of :math:`{\mu^*}` or its rank with the number of trajectories. Parameters @@ -1273,7 +1273,7 @@ def plot_morris_convergence(result_dct, metric, parameters=(), ylabel = f'Rank for {metric.name.lower()}' loc = 'lower left' else: - ylabel = f'$\mu^*$ for {metric.name.lower()}' + ylabel = rf'$\mu^*$ for {metric.name.lower()}' loc = 'best' palette = sns.color_palette('deep', n_colors=len(param_names)) @@ -1357,7 +1357,7 @@ def _plot_heatmap(hmap_df, ax=None, annot=False, diagonal='', sts1_df=None, def plot_fast_results(result_dct, metric, parameters=(), ax=None, error_bar=True, file='', close_fig=True): - ''' + r''' Visualize the results from FAST or RBD-FAST analysis as a bar plot. Parameters @@ -1410,7 +1410,7 @@ def plot_sobol_results(result_dct, metric, ax=None, parameters=(), kind='all', annotate_heatmap=False, plot_in_diagonal='', error_bar=True, file='', close_fig=True): - ''' + r''' Visualize the results from Sobol analysis as a bar plot and/or heat map. Total (:math:`S_{Ti}`) and main (:math:`S_{1i}`) effects can be drawn in the bar plot or diagonal of the heat map; From ccb7314a61f7c5813fc9b59a594d444edf87480e Mon Sep 17 00:00:00 2001 From: Yalin Li Date: Tue, 10 Dec 2024 14:27:33 -0500 Subject: [PATCH 27/32] temporary fix for indexer phase (thermosteam recent updates) --- qsdsan/_waste_stream.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/qsdsan/_waste_stream.py b/qsdsan/_waste_stream.py index 37c5843d..4700c430 100644 --- a/qsdsan/_waste_stream.py +++ b/qsdsan/_waste_stream.py @@ -154,7 +154,11 @@ def __init__(self, dct, F_vol, MW, phase, phase_container): def output(self, index, value): '''Concentration flows, in mg/L (g/m3).''' f_mass = value * self.MW[index] - phase = self.phase or self.phase_container.phase + if self.phase: + phase = self.phase + else: + try: phase = self.phase_container._phase + except: phase = self.phase_container if phase != 'l': raise AttributeError('Concentration only valid for liquid phase.') V_sum = self.F_vol @@ -186,7 +190,7 @@ def by_conc(self, TP): check_data=False, ) return conc -indexer.ChemicalMolarFlowIndexer.by_conc = by_conc +ChemicalMolarFlowIndexer.by_conc = by_conc del by_conc From 3f40a30e999475e0e3a1b123ebd79563d37c5000 Mon Sep 17 00:00:00 2001 From: Yalin Li Date: Tue, 10 Dec 2024 14:31:10 -0500 Subject: [PATCH 28/32] add new attribute with biosteam update --- qsdsan/_sanunit.py | 1 + 1 file changed, 1 insertion(+) diff --git a/qsdsan/_sanunit.py b/qsdsan/_sanunit.py index 51cac1a7..4788b965 100644 --- a/qsdsan/_sanunit.py +++ b/qsdsan/_sanunit.py @@ -217,6 +217,7 @@ def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='WasteStream self._assert_compatible_property_package() self._utility_cost = None + self._recycle_system = None ##### qsdsan-specific ##### for i in (*construction, *transportation, *equipment): From d2edfe64656a29e3a83acd374588f7b0d71a9ad9 Mon Sep 17 00:00:00 2001 From: Yalin Li Date: Tue, 10 Dec 2024 15:39:13 -0500 Subject: [PATCH 29/32] update version --- docs/source/conf.py | 2 +- setup.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index cc735083..20cc44a7 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -31,7 +31,7 @@ # built documents. # # The short X.Y version. -version = '1.4.0' +version = '1.4.1' # The full version, including alpha/beta/rc tags. release = version diff --git a/setup.py b/setup.py index 2059114b..995ef6d5 100644 --- a/setup.py +++ b/setup.py @@ -18,7 +18,7 @@ setup( name='qsdsan', packages=['qsdsan'], - version='1.4.0', + version='1.4.1', license='UIUC', author='Quantitative Sustainable Design Group', author_email='quantitative.sustainable.design@gmail.com', From 557cade22cafb91191cb41dab858af5e5b12d973 Mon Sep 17 00:00:00 2001 From: Yalin Li Date: Tue, 10 Dec 2024 17:06:37 -0500 Subject: [PATCH 30/32] resolve conflict --- requirements.txt | 3 +++ 1 file changed, 3 insertions(+) diff --git a/requirements.txt b/requirements.txt index c9e85842..8fa8076c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,7 +6,10 @@ SALib>=1.4.5 seaborn sympy>=1.8 matplotlib<=3.6.0 +<<<<<<< Updated upstream # numpy<2.0.0 +======= +>>>>>>> Stashed changes # Specifically for docs sphinx From 0c70f425e05737700d5179b9b5a9fa64f9a9eb88 Mon Sep 17 00:00:00 2001 From: Yalin Li Date: Tue, 10 Dec 2024 17:07:45 -0500 Subject: [PATCH 31/32] remove previous `matplotlib` version requirement --- requirements.txt | 4 ---- 1 file changed, 4 deletions(-) diff --git a/requirements.txt b/requirements.txt index 8fa8076c..37a3c2d7 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,10 +6,6 @@ SALib>=1.4.5 seaborn sympy>=1.8 matplotlib<=3.6.0 -<<<<<<< Updated upstream -# numpy<2.0.0 -======= ->>>>>>> Stashed changes # Specifically for docs sphinx From e45b7b43fbb4f56ff18158ad05972b1cc7ecbfd7 Mon Sep 17 00:00:00 2001 From: Yalin Li Date: Tue, 10 Dec 2024 17:08:16 -0500 Subject: [PATCH 32/32] remove previous `matplotlib` version requirement --- requirements.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 37a3c2d7..22b2a295 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,7 +5,6 @@ git+https://github.com/QSD-Group/EXPOsan.git SALib>=1.4.5 seaborn sympy>=1.8 -matplotlib<=3.6.0 # Specifically for docs sphinx