From 9afc0339ce19643020a52a51a8fb62a14910f683 Mon Sep 17 00:00:00 2001 From: Yalin Date: Wed, 16 Oct 2024 07:56:27 -0400 Subject: [PATCH 01/12] update doc version --- docs/source/conf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index 7ea627c7..cc735083 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -31,7 +31,7 @@ # built documents. # # The short X.Y version. -version = '1.3.1' +version = '1.4.0' # The full version, including alpha/beta/rc tags. release = version From a2801feda5b3eab634f3868ecc9bfa6cd01e5067 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Mon, 21 Oct 2024 10:58:59 -0700 Subject: [PATCH 02/12] add tests for junctions --- tests/test_junctions.py | 497 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 497 insertions(+) create mode 100644 tests/test_junctions.py diff --git a/tests/test_junctions.py b/tests/test_junctions.py new file mode 100644 index 00000000..3edad88b --- /dev/null +++ b/tests/test_junctions.py @@ -0,0 +1,497 @@ +# -*- coding: utf-8 -*- +''' +EXPOsan: Exposition of sanitation and resource recovery systems + +This module is developed by: + + Joy Zhang + + Yalin Li + +This module is under the University of Illinois/NCSA Open Source License. +Please refer to https://github.com/QSD-Group/EXPOsan/blob/main/LICENSE.txt +for license details. + +Reference: +.. [1] Alex, J.; Benedetti, L.; Copp, J. B.; Gernaey, K. V.; Jeppsson, U.; + Nopens, I.; Pons, M. N.; Rosen, C.; Steyer, J. P.; Vanrolleghem, P. A. + Benchmark Simulation Model No. 2 (BSM2). + http://iwa-mia.org/wp-content/uploads/2022/09/TR3_BSM_TG_Tech_Report_no_3_BSM2_General_Description.pdf. +.. [2] Flores-Alsina, X., Solon, K., Kazadi Mbamba, C., Tait, S., Gernaey, K. v., + Jeppsson, U., & Batstone, D. J. (2016). Modelling phosphorus (P), + sulfur (S) and iron (Fe) interactions for dynamic simulations of anaerobic + digestion processes. Water Research, 95, 370–382. + https://doi.org/10.1016/J.WATRES.2016.03.012 + +''' +#%% + +def test_adm1_junctions(): + + import qsdsan as qs, numpy as np + from numpy.testing import assert_allclose as ac + from qsdsan import ( + processes as pc, + sanunits as su, + WasteStream, + ) + from qsdsan.utils import ospath, load_data + from exposan.bsm2 import data_path + + matlab_preAD_adm = { + 'S_su': 0.0, # monosacharides (kg COD/m3) + 'S_aa': 0.04388, # amino acids (kg COD/m3) + 'S_fa': 0.0, # long chain fatty acids (LCFA) (kg COD/m3) + 'S_va': 0.0, # total valerate (kg COD/m3) + 'S_bu': 0.0, # total butyrate (kg COD/m3) + 'S_pro': 0.0, # total propionate (kg COD/m3) + 'S_ac': 0.0, # total acetate (kg COD/m3) + 'S_h2': 0.0, # hydrogen gas (kg COD/m3) + 'S_ch4': 0.0, # methane gas (kg COD/m3) + 'S_IC': 0.0079326*12, # inorganic carbon (kmole C/m3 -> kg C/m3) 0.0951912 + 'S_IN': 0.0019721*14, # inorganic nitrogen (kmole N/m3 -> kg N/m3) 0.0276094 + 'S_I': 0.028067, # soluble inerts (kg COD/m3) + 'X_c': 0.0, # composites (kg COD/m3) + 'X_ch': 3.7236, # carbohydrates (kg COD/m3) + 'X_pr': 15.9235, # proteins (kg COD/m3) + 'X_li': 8.047, # lipids (kg COD/m3) + 'X_su': 0.0, # sugar degraders (kg COD/m3) + 'X_aa': 0.0, # amino acid degraders (kg COD/m3) + 'X_fa': 0.0, # LCFA degraders (kg COD/m3) + 'X_c4': 0.0, # valerate and butyrate degraders (kg COD/m3) + 'X_pro': 0.0, # propionate degraders (kg COD/m3) + 'X_ac': 0.0, # acetate degraders (kg COD/m3) + 'X_h2': 0.0, # hydrogen degraders (kg COD/m3) + 'X_I': 17.0106, # particulate inerts (kg COD/m3) + 'S_cat': 0.0, # cations (base) (kmole/m3) + 'S_an': 0.0052101, # anions (acid) (kmole/m3) + # 'Q': 178.4674, # Flow rate (m3/d) + } + + matlab_postAD_adm = { + 'S_su': 0.012394, + 'S_aa': 0.0055432, + 'S_fa': 0.10741, + 'S_va': 0.012333, + 'S_bu': 0.014003, + 'S_pro': 0.017584, + 'S_ac': 0.089315, + 'S_h2': 2.5055e-07, + 'S_ch4': 0.05549, + 'S_IC': 0.095149*12, + 'S_IN': 1.3226, + 'S_I': 0.13087, + 'X_c': 0.10792, + 'X_ch': 0.020517, + 'X_pr': 0.08422, + 'X_li': 0.043629, + 'X_su': 0.31222, + 'X_aa': 0.93167, + 'X_fa': 0.33839, + 'X_c4': 0.33577, + 'X_pro': 0.10112, + 'X_ac': 0.67724, + 'X_h2': 0.28484, + 'X_I': 17.2162, + 'S_cat': 0., #-4.0789e-34, + 'S_an': 0.0052101 + } + + matlab_postAD_asm = { + 'S_I': 130.867, # soluble inert organic matter, mg COD/l + 'S_S': 258.5789, # readily biodegradable substrate, mg COD/l + 'X_I': 17216.2434, # particulate inert organic matter, mg COD/l + 'X_S': 2611.4843, # slowly biodegradable substrate, mg COD/l + 'X_BH': 0.0, # active heterotrophic biomass, mg COD/l + 'X_BA': 0.0, # active autotrophic biomass, mg COD/l + 'X_P': 626.0652, # particulate products arising from biomass decay, mg COD/l + 'S_O': 0.0, # dissolved O2, mg -COD/l + 'S_NO': 0.0, # nitrate and nitrite nitrogen, mg N/L + 'S_NH': 1442.7882, # ammonium, mg N/L + 'S_ND': 0.54323, # soluble biodegradable organic nitrogen + 'X_ND': 100.8668, # particulate biodegradable organic nitrogen, mg N/l + 'S_ALK': 97.8459*12, # alkalinity, assumed to be HCO3-, 97.8459, mol HCO3/m3 -> g C/m3 + 'S_N2': 0.0, # dissolved O2 + # 'Q': 178.4674, # Flow rate, m3/d + } + + + adm1init = load_data(ospath.join(data_path, 'adm1init.csv'), index_col=0).to_dict('index') + asm1_default_parameters = dict( + mu_H = 4.0, + K_S = 10.0, + K_OH = 0.2, + K_NO = 0.5, + b_H = 0.3, + mu_A = 0.5, + K_NH = 1.0, + K_OA = 0.4, + b_A = 0.05, + eta_g = 0.8, + k_a = 0.05, + k_h = 3.0, + K_X = 0.1, + eta_h = 0.8, + Y_H = 0.67, + Y_A = 0.24, + f_P = 0.08, + i_XB = 0.08, + i_XP = 0.06, + fr_SS_COD = 0.75 + ) + + T = 273.15 + 35 + cmps_asm1 = pc.create_asm1_cmps() + asm1 = pc.ASM1(components=cmps_asm1, **asm1_default_parameters) + preAD_asm = WasteStream('preAD_asm', T=T) + preAD_asm.set_flow_by_concentration( + flow_tot=178.4674, + concentrations=dict( + S_I = 28.0665, + S_S = 48.9526, + X_I = 10361.7101, + X_S = 20375.0176, + X_BH = 10210.0698, + X_BA = 553.2808, + X_P = 3204.6601, + S_O = 0.25225, + S_NO = 1.6871, + S_NH = 28.9098, + S_ND = 4.6834, + X_ND = 906.0933, + S_ALK = 7.1549*12 + ), + units=('m3/d', 'mg/L') + ) + thermo_asm1 = qs.get_thermo() + cmps_adm1 = pc.create_adm1_cmps() + adm1 = pc.ADM1() + cmps_adm1.X_I.i_N = cmps_asm1.X_I.i_N # slight difference + cmps_adm1.refresh_constants() + thermo_adm1 = qs.get_thermo() + + J1 = su.ASMtoADM('J1', upstream=preAD_asm, downstream='preAD_adm', + thermo=thermo_adm1, isdynamic=True, adm1_model=adm1,#) + T=T, pH=7.2631) + AD1 = su.AnaerobicCSTR('AD1', ins=J1-0, outs=('biogas', 'postAD_adm'), + isdynamic=True, V_liq=3400, V_gas=300, T=T, + model=adm1,) + AD1.set_init_conc(**adm1init['AD1']) + # Switch back to ASM1 components + J2 = su.ADMtoASM('J2', upstream=AD1-1, downstream='postAD_asm', + thermo=thermo_asm1, isdynamic=True, adm1_model=adm1) + J2.bio_to_xs = 0.79 + qs.set_thermo(thermo_asm1) + + sys = qs.System(path=(J1, AD1, J2)) + sys.simulate(state_reset_hook='reset_cache', t_span=(0, 200), method='BDF') + fs = sys.flowsheet.stream + + for ws in sys.streams: + ws.state[ws.state < 2.2e-16] = 0 + + ac(cmps_adm1.kwarray(matlab_preAD_adm)[:-1]*1e3, fs.preAD_adm.state[:-2], rtol=1e-4) + ac(cmps_adm1.kwarray(matlab_postAD_adm)[:-1]*1e3, fs.postAD_adm.state[:-2], rtol=1e-2) + ac(cmps_asm1.kwarray(matlab_postAD_asm)[:-1], fs.postAD_asm.state[:-2], rtol=1e-3) + + h2 = cmps_adm1.S_h2 + ch4 = cmps_adm1.S_ch4 + co2 = cmps_adm1.S_IC + assert np.isclose(AD1.state['S_h2_gas'] * h2.chem_MW / h2.i_mass, 1.1032e-5, rtol=1e-3) + assert np.isclose(AD1.state['S_ch4_gas'] * ch4.chem_MW / ch4.i_mass, 1.6535, rtol=1e-2) + assert np.isclose(AD1.state['S_IC_gas'], 0.01354, rtol=1e-2) + assert np.isclose(AD1.outs[1].pH, 7.2631, rtol=1e-3) + + assert np.isclose(fs.biogas.imass['S_h2']*24 * h2.i_mass, 0.0035541, rtol=1e-2) + assert np.isclose(fs.biogas.imass['S_ch4']*24 * ch4.i_mass, 1065.3523, rtol=1e-2) + assert np.isclose(fs.biogas.imass['S_IC']*24 * co2.i_mass, 1535.4118, rtol=1e-2) + + sys.flowsheet.clear() + +#%% +def test_adm1p_junctions(): + import numpy as np + from numpy.testing import assert_allclose as ac + from chemicals.elements import molecular_weight as get_mw + from qsdsan import sanunits as su, processes as pc, WasteStream, System, get_thermo + # from qsdsan.utils import load_data, ospath, time_printer + # from exposan.bsm2 import data_path + + Q = 190 # influent flowrate [m3/d] + HRT = 20 + V_liq = Q*HRT + V_gas = 0.088*V_liq + Temp = 273.15+35 # temperature [K] + C_mw = get_mw({'C':1}) + N_mw = get_mw({'N':1}) + P_mw = get_mw({'P':1}) + struv_mw = get_mw(dict(Mg=1, N=1, H=4, P=1, O=4)) + # adm1init = load_data(ospath.join(data_path, 'adm1init.csv'), index_col=0).to_dict('index') + + # Table 1.1 [mg/L], Flores-Alsina et al., 2016. Appendix + inf_asm2d = dict( + S_O2=0, + S_F=26.44, + S_A=17.66, + S_I=27.23, + S_NH4=18.58, + S_N2=5.07, + S_NO3=0.02, + S_PO4=4.69, + S_IC=78.99, + X_I=10964.41, + X_S=19084.76, + X_H=9479.39, + X_PAO=3862.20, + X_PP=450.87, + X_PHA=24.64, + X_AUT=333.79, + S_K=19.79, + S_Mg=189.87, + S_Na=70, + S_Cl=1035, + S_Ca=300, + ) + + # Table 1.3 [kg/m3] + inf_adm1p = dict( + S_su=0.018, + S_aa=0.008, + S_ac=0.018, + S_IC=0.021*C_mw, + S_IN=0.036*N_mw, + S_IP=0.006*P_mw, + S_I=0.027, + X_ch=8.020, + X_pr=8.481, + X_li=11.416, + X_I=11.946, + X_PHA=0.025, + X_PP=0.015*P_mw, + X_PAO=3.862, + S_K=0.001*39, + S_Mg=0.008*24.3, + S_Ca=0.007*40, + S_Na=0.003*23, + S_Cl=0.029*35.5, + # S_N2=0.0004*14 + ) + + # [kmol/m3] + _inf_adm1p = dict( + S_IC=0.021, + S_IN=0.036, + S_IP=0.006, + X_PP=0.015, + S_K=0.001, + S_Mg=0.008, + S_Ca=0.007, + S_Na=0.003, + S_Cl=0.029, + ) + + # Table 1.4 [kg/m3] + out_adm1p = dict( + S_su=0.013, + S_aa=0.006, + S_fa=0.116, + S_va=0.012, + S_bu=0.016, + S_pro=0.019, + S_ac=0.055, + S_h2=2.65e-7, + S_ch4=0.052, + S_IC=0.059*C_mw, + S_IN=0.080*N_mw, + S_IP=0.007*P_mw, + S_I=0.027, + X_ch=1.441, + X_pr=1.513, + X_li=2.025, + X_I=12.345, + X_PHA=0.252, + X_PP=8.05e-6*P_mw, + # X_biomass=3.600, + X_su=3.600, + S_K=0.005*39, + S_Mg=0.001*24.3, + S_Ca=0.001*40, + X_ACP=0.002*310.176722, + X_struv=0.011*245.406502, + S_Na=0.003*23, + S_Cl=0.029*35.5, + # S_N2=0.0004*14 + ) + + # _out_adm1p = dict( + # S_IC=0.059, + # S_IN=0.080, + # S_IP=0.007, + # X_PP=8.05e-6, + # S_K=0.005, + # S_Mg=0.001, + # S_Ca=0.001, + # X_ACP=0.002, + # X_struv=0.011, + # S_Na=0.003, + # S_Cl=0.029, + # ) + + # Table 1.5 [mg/L] + out_asm2d = dict( + S_NH4=1291.68, + S_PO4=298.09, + S_F=134.43, + S_A=353.82, + S_I=27.23, + S_IC=885.27, + S_K=208.84, + S_Mg=28.29, + X_I=12704.93, + X_S=8218.94, + S_Na=70, + S_Cl=1035, + S_Ca=20.45, + X_ACP=722.17, + X_struv=1578.52*245.406502/struv_mw + ) + + # [mmol/L] + _out_asm2d = dict( + S_NH4=1291.68/N_mw, + S_PO4=298.09/P_mw, + S_IC=885.27/C_mw, + S_K=208.84/39, + S_Mg=28.29/24.3, + S_Na=70/23, + S_Cl=1035/35.5, + S_Ca=20.45/40, + X_ACP=722.17/310.176722, + X_struv=1578.52/struv_mw + ) + + default_init_conds = { + 'S_su': 0.014*1e3, + 'S_aa': 0.0062*1e3, + 'S_fa': 0.126*1e3, + 'S_va': 0.0129*1e3, + 'S_bu': 0.0168*1e3, + 'S_pro': 0.0204*1e3, + 'S_ac': 0.0588*1e3, + 'S_h2': 2.8309e-7*1e3, + 'S_ch4': 0.0544*1e3, + 'S_IC': 0.089*12*1e3, + 'S_IN': 0.0663*14*1e3, + 'S_IP': 0.028*31*1e3, + 'S_I': 0.1309*1e3, + 'X_ch': 1.302*1e3, + 'X_pr': 1.3613*1e3, + 'X_li': 1.8127*1e3, + 'X_su': 0.5146*1e3, + 'X_aa': 0.4017*1e3, + 'X_fa': 0.3749*1e3, + 'X_c4': 0.1596*1e3, + 'X_pro': 0.0896*1e3, + 'X_ac': 0.5006*1e3, + 'X_h2': 0.258*1e3, + 'X_I': 12.9232*1e3, + 'X_PHA': 0.6697*1e3, + 'X_PAO': 0.9154*1e3, + 'S_K': 0.0129*1e3, + 'S_Mg': 0.0001*1e3, + 'S_Ca': 2e-4*1e3, + 'X_struv':0.0161*1e3, + 'X_ACP': 9e-4*1e3, + 'X_FePO4': 0.001*1e3, + 'S_Na': 0.061*1e3, + 'S_Cl': 0.0126*1e3 + } + + cmps_asm = pc.create_masm2d_cmps() + inf_asm = WasteStream('inf_asm', T=Temp) + inf_asm.set_flow_by_concentration( + flow_tot=Q, + concentrations=inf_asm2d, + units=('m3/d', 'mg/L') + ) + alt_eff_asm = WasteStream('alt_eff_asm', T=Temp) + alt_eff_asm.set_flow_by_concentration( + flow_tot=Q, + concentrations=out_asm2d, + units=('m3/d', 'mg/L') + ) + asm = pc.mASM2d() + thermo_asm = get_thermo() + cmps_adm = pc.create_adm1p_cmps() + alt_inf_adm = WasteStream('alt_inf_adm', T=Temp) + alt_inf_adm.set_flow_by_concentration( + flow_tot=Q, + concentrations=inf_adm1p, + units=('m3/d', 'kg/m3') + ) + alt_eff_adm = WasteStream('alt_eff_adm', T=Temp) + alt_eff_adm.set_flow_by_concentration( + flow_tot=Q, + concentrations=out_adm1p, + units=('m3/d', 'kg/m3') + ) + adm = pc.ADM1p( + f_bu_su=0.1328, f_pro_su=0.2691, f_ac_su=0.4076, + q_ch_hyd=0.3, q_pr_hyd=0.3, q_li_hyd=0.3, + ) + thermo_adm = get_thermo() + + J1 = su.mASM2dtoADM1p( + 'J1', upstream=inf_asm, downstream='inf_adm', + thermo=thermo_adm, isdynamic=True, + adm1_model=adm, asm2d_model=asm + ) + J1.xs_to_li = 0.6 + AD = su.AnaerobicCSTR( + 'AD', + ins=alt_inf_adm, + # ins=J1-0, + outs=('biogas', 'eff_adm'), isdynamic=True, + V_liq=V_liq, V_gas=V_gas, T=Temp, model=adm + ) + AD.algebraic_h2 = False + AD.set_init_conc(**default_init_conds) + J2 = su.ADM1ptomASM2d( + 'J2', + upstream=alt_eff_adm, + # upstream=AD-1, + downstream='eff_asm', thermo=thermo_asm, isdynamic=True, + adm1_model=adm, asm2d_model=asm + ) + + sys = System(path=(J1, AD, J2)) + sys.simulate(state_reset_hook='reset_cache', t_span=(0, 200), method='BDF') + s = sys.flowsheet.stream + + ########## mASM2d to ADM1p ########### + mass2mol = cmps_adm.i_mass / cmps_adm.chem_MW + idx = cmps_adm.indices(_inf_adm1p.keys()) + _molar = np.round(s.inf_adm.conc[idx] * mass2mol[idx] * 1e-3, 3) + ac(_molar, np.array([v for v in _inf_adm1p.values()])) + ac(np.delete(s.inf_adm.conc, idx)[:-1], # exclude water + np.delete(s.alt_inf_adm.conc, idx)[:-1], + atol=1.0) + + ########## !!! ADM1p skip for now ########## + + ######### ADM1p to mASM2d ########### + mass2mol = cmps_asm.i_mass / cmps_asm.chem_MW + idx = cmps_asm.indices(_out_asm2d.keys()) + _molar = s.eff_asm.conc[idx] * mass2mol[idx] + ac(_molar, np.array([v for v in _out_asm2d.values()]), atol=1.0) + ac(np.delete(s.eff_asm.conc, idx)[:-1], # exclude water + np.delete(s.alt_eff_asm.conc, idx)[:-1], + atol=1.0) + + sys.flowsheet.clear() + +#%% + +if __name__ == '__main__': + test_adm1_junctions() + test_adm1p_junctions() From 564a47886eb87b6d4a5a6c4e2529301185d8ae62 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Tue, 22 Oct 2024 10:51:13 -0700 Subject: [PATCH 03/12] add doctest for `DiffusedAeration` --- qsdsan/processes/_aeration.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/qsdsan/processes/_aeration.py b/qsdsan/processes/_aeration.py index d5a0ffec..9da6d678 100644 --- a/qsdsan/processes/_aeration.py +++ b/qsdsan/processes/_aeration.py @@ -80,6 +80,25 @@ class DiffusedAeration(Process): [parameters] KLa: 240 DOsat: 8 [dynamic parameters] + + >>> aer2 = pc.DiffusedAeration('aer2', 'S_O', KLa_20=100, V=1000, d_submergence=3.7) + >>> aer2.show() + Process: aer2 + [stoichiometry] S_O: 1 + [reference] S_O + [rate equation] KLa*(DOsat - S_O) + [parameters] KLa: 60 + DOsat: 9.87 + [dynamic parameters] + + >>> aer2.Q_air # doctest: +ELLIPSIS + 12470.65... + >>> round(aer2.SOTR / 1000) + 1039 + + >>> aer3 = pc.DiffusedAeration('aer3', 'S_O', Q_air=7600, V=1000, d_submergence=3.7) + >>> aer3.kLa # doctest: +ELLIPSIS + 36.56... """ @@ -357,5 +376,7 @@ def DOsat(self): return self._DOsat @DOsat.setter def DOsat(self, DOsat): + if DOsat is not None: + self._DOsat_s20 = DOsat / (self.tau * self._beta * self.Omega * self.delta) self._DOsat = DOsat or self._calc_DOsat() self.set_parameters(DOsat = self._DOsat) \ No newline at end of file From b8b55aa14ed392dfdb84c58087e7458a77f49c6e Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Tue, 22 Oct 2024 14:14:23 -0700 Subject: [PATCH 04/12] add test for `mASM2d` --- qsdsan/processes/_asm2d.py | 46 +++++++++++++++++++++++++++++++++++--- 1 file changed, 43 insertions(+), 3 deletions(-) diff --git a/qsdsan/processes/_asm2d.py b/qsdsan/processes/_asm2d.py index d12e90c5..c2a99b7e 100644 --- a/qsdsan/processes/_asm2d.py +++ b/qsdsan/processes/_asm2d.py @@ -73,8 +73,6 @@ def create_asm2d_cmps(set_thermo=True): return cmps_asm2d -# create_asm2d_cmps() - def create_masm2d_cmps(set_thermo=True): c2d = create_asm2d_cmps(False) ion_kwargs = dict(particle_size='Soluble', @@ -433,7 +431,7 @@ class ASM2d(CompiledProcesses): Examples -------- - >>> from qsdsan import processes as pc, set_thermo + >>> from qsdsan import processes as pc >>> cmps = pc.create_asm2d_cmps() >>> asm2d = pc.ASM2d() >>> asm2d.show() @@ -714,6 +712,8 @@ class mASM2d(CompiledProcesses): electron_acceptor_dependent_decay : bool, optional Whether biomass decay kinetics is dependent on concentrations of electron acceptors. The default is True. + pH_ctrl : float or None, optional + Whether to fix pH at a specific value or solve for pH (`None`). The default is 7.0. k_h : float, optional Hydrolysis rate constant, in [d^(-1)]. The default is 3.0. eta_NO3_Hl : float, optional @@ -766,6 +766,46 @@ class mASM2d(CompiledProcesses): >>> asm.show() mASM2d([aero_hydrolysis, anox_hydrolysis, anae_hydrolysis, hetero_growth_S_F, hetero_growth_S_A, denitri_S_F, denitri_S_A, ferment, hetero_lysis, storage_PHA, aero_storage_PP, anox_storage_PP, PAO_aero_growth_PHA, PAO_anox_growth, PAO_lysis, PP_lysis, PHA_lysis, auto_aero_growth, auto_lysis, CaCO3_precipitation_dissolution, struvite_precipitation_dissolution, newberyite_precipitation_dissolution, ACP_precipitation_dissolution, MgCO3_precipitation_dissolution, AlPO4_precipitation_dissolution, FePO4_precipitation_dissolution]) + >>> # Calculate process rate given state variable values and fixed pH. + >>> import numpy as np + >>> state_arr = np.ones(len(cmps)) + >>> rhos = asm.rate_function(state_arr) # reaction rate for each process + >>> for i,j in zip(asm.IDs, rhos): + ... print(f'{i}{(40-len(i))*" "}{j:.3g}') + aero_hydrolysis 2.27 + anox_hydrolysis 0.182 + anae_hydrolysis 0.0606 + hetero_growth_S_F 0.471 + hetero_growth_S_A 0.471 + denitri_S_F 0.0503 + denitri_S_A 0.0503 + ferment 0.0333 + hetero_lysis 0.356 + storage_PHA 0.594 + aero_storage_PP 1.06 + anox_storage_PP 0.0851 + PAO_aero_growth_PHA 0.778 + PAO_anox_growth 0.0622 + PAO_lysis 0.174 + PP_lysis 0.174 + PHA_lysis 0.174 + auto_aero_growth 0.33 + auto_lysis 0.111 + CaCO3_precipitation_dissolution 0 + struvite_precipitation_dissolution 0 + newberyite_precipitation_dissolution 0 + ACP_precipitation_dissolution 0 + MgCO3_precipitation_dissolution 0 + AlPO4_precipitation_dissolution 1.82e-11 + FePO4_precipitation_dissolution 1.82e-11 + + >>> # Estimate pH given state variable values. + >>> Ka = asm.rate_function.params['Ka'] + >>> unit_conversion = asm.rate_function.params['mass2mol'] + >>> h_ion = asm.solve_pH(state_arr, Ka, unit_conversion) + >>> pH = -np.log10(h_ion) + >>> print(f'{pH:.2f}') + 8.40 References ---------- From 7f566bb120acfa26353304cf43552485f65d2c29 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Tue, 22 Oct 2024 15:30:23 -0700 Subject: [PATCH 05/12] added doctest for `ADM1_p_extension` --- qsdsan/processes/_adm1_p_extension.py | 37 ++++++++++++++++++++++++++- 1 file changed, 36 insertions(+), 1 deletion(-) diff --git a/qsdsan/processes/_adm1_p_extension.py b/qsdsan/processes/_adm1_p_extension.py index 1c1a81cb..ab342969 100644 --- a/qsdsan/processes/_adm1_p_extension.py +++ b/qsdsan/processes/_adm1_p_extension.py @@ -239,7 +239,6 @@ def _rhos_adm1_p_extension(state_arr, params, h=None): rhos[9] *= Inh3 rhos[-3:] = kLa * (biogas_S - KH * biogas_p) - # print(rhos) return rhos def dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, V_liq, S_h2_in): @@ -355,6 +354,42 @@ class ADM1_p_extension(ADM1): >>> adm1_p.show() ADM1_p_extension([hydrolysis_carbs, hydrolysis_proteins, hydrolysis_lipids, uptake_sugars, uptake_amino_acids, uptake_LCFA, uptake_valerate, uptake_butyrate, uptake_propionate, uptake_acetate, uptake_h2, decay_Xsu, decay_Xaa, decay_Xfa, decay_Xc4, decay_Xpro, decay_Xac, decay_Xh2, storage_Sva_in_XPHA, storage_Sbu_in_XPHA, storage_Spro_in_XPHA, storage_Sac_in_XPHA, lysis_XPAO, lysis_XPP, lysis_XPHA, h2_transfer, ch4_transfer, IC_transfer]) + >>> import numpy as np + >>> state_arr = np.ones(cmps.size + len(adm1_p._biogas_IDs) + 2) # liquid-phase concentrations, gas-phase concentrations, liquid flowrate, and temperature + >>> state_arr[-1] = 273.15+35 # Temperature + >>> rhos = adm1_p.rate_function(state_arr) # reaction rate for each process + >>> for i,j in zip(adm1_p.IDs, rhos): + ... print(f'{i}{(40-len(i))*" "}{j:.3g}') + hydrolysis_carbs 10 + hydrolysis_proteins 10 + hydrolysis_lipids 10 + uptake_sugars 20 + uptake_amino_acids 38.4 + uptake_LCFA 2.14e-05 + uptake_valerate 8.32e-05 + uptake_butyrate 8.32e-05 + uptake_propionate 4.13e-05 + uptake_acetate 1.93 + uptake_h2 34.9 + decay_Xsu 0.02 + decay_Xaa 0.02 + decay_Xfa 0.02 + decay_Xc4 0.02 + decay_Xpro 0.02 + decay_Xac 0.02 + decay_Xh2 0.02 + storage_Sva_in_XPHA 0.747 + storage_Sbu_in_XPHA 0.747 + storage_Spro_in_XPHA 0.747 + storage_Sac_in_XPHA 0.747 + lysis_XPAO 0.2 + lysis_XPP 0.2 + lysis_XPHA 0.2 + h2_transfer 139 + ch4_transfer -181 + IC_transfer -1.66e+03 + + References ---------- [1] Batstone, D. J.; Keller, J.; Angelidaki, I.; Kalyuzhnyi, S. V; From 69140a4e81f523000ebf0197541f9094da1a7cd8 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Wed, 23 Oct 2024 09:52:37 -0700 Subject: [PATCH 06/12] add doctest for `BatchExperiment` --- qsdsan/processes/_adm1_p_extension.py | 9 - .../sanunits/_suspended_growth_bioreactor.py | 513 +++++++++--------- 2 files changed, 269 insertions(+), 253 deletions(-) diff --git a/qsdsan/processes/_adm1_p_extension.py b/qsdsan/processes/_adm1_p_extension.py index ab342969..aa4d1818 100644 --- a/qsdsan/processes/_adm1_p_extension.py +++ b/qsdsan/processes/_adm1_p_extension.py @@ -94,15 +94,6 @@ def acid_base_rxn(h_ion, weak_acids_tot, Kas): nh3, hpo4, hco3, ac, pro, bu, va = Kas[1:] * weak_acids_tot[4:] / (Kas[1:] + h_ion) return S_cat + S_K + 2*S_Mg + h_ion + (S_IN - nh3) - S_an - oh_ion - hco3 - ac - pro - bu - va - 2*hpo4 - (S_IP - hpo4) -# The function 'fprime_abr' is not used in the code -def fprime_abr(h_ion, weak_acids_tot, Kas): - # S_cat, S_K, S_Mg, S_an, S_IN, S_IP = weak_acids_tot[:6] - Kw = Kas[0] - doh_ion = - Kw / h_ion ** 2 - dnh3, dhpo4, dhco3, dac, dpro, dbu, dva = - Kas[1:] * weak_acids_tot[4:] / (Kas[1:] + h_ion)**2 - return 1 + (-dnh3) - doh_ion - dhco3 - dac - dpro - dbu - dva - dhpo4 - - rhos = np.zeros(28) # 28 kinetic processes (25 as defined in modified ADM1 + 3 for gases) Cs = np.empty(25) # 25 processes as defined in modified ADM1 diff --git a/qsdsan/sanunits/_suspended_growth_bioreactor.py b/qsdsan/sanunits/_suspended_growth_bioreactor.py index 5c305a59..0196daf8 100644 --- a/qsdsan/sanunits/_suspended_growth_bioreactor.py +++ b/qsdsan/sanunits/_suspended_growth_bioreactor.py @@ -10,20 +10,19 @@ for license details. ''' -from .. import SanUnit, WasteStream, Process, Processes, CompiledProcesses -from ._clarifier import _settling_flux +from .. import SanUnit, WasteStream, Process, CompiledProcesses +# from ._clarifier import _settling_flux from ..sanunits import dydt_cstr -from sympy import symbols, lambdify, Matrix -from scipy.integrate import solve_ivp +# from scipy.integrate import solve_ivp from warnings import warn -from math import floor, ceil +# from math import floor, ceil import numpy as np import pandas as pd -from numba import njit +# from numba import njit __all__ = ('CSTR', 'BatchExperiment', - 'SBR', + # 'SBR', 'PFR', ) @@ -437,7 +436,42 @@ def _cost(self): #%% class BatchExperiment(SanUnit): + ''' + A batch reactor in experimental settings. + + Parameters + ---------- + model : :class:`CompiledProcesses`, optional + Process model that describes the dynamics of state variables. + The `state` of the batch reactor is entirely determined by the + stoichiometry and rate function in this model. + + Examples + -------- + >>> import qsdsan.sanunits as su, qsdsan.processes as pc + >>> cmps = pc.create_asm1_cmps() + >>> asm1 = pc.ASM1() + >>> BE = su.BatchExperiment('BE', model=asm1) + >>> BE.set_init_conc(S_S=20, X_BH=500, S_O=8, S_ND=3, S_ALK=84) + >>> BE.simulate(t_span=(0,10), method='BDF') + >>> BE.state + {'S_I': 0.0, + 'S_S': 0.9259111676784925, + 'X_I': 0.0, + 'X_S': 446.163046741642, + 'X_BH': 25.664835161702214, + 'X_BA': 0.0, + 'X_P': 39.246207146417575, + 'S_O': -1.7987340043125508e-08, + 'S_NO': -4.085240753452882e-20, + 'S_NH': 2.1230597299994245, + 'S_ND': 1.6652546123923254e-06, + 'X_ND': 36.46897947131222, + 'S_ALK': 85.82051689686057, + 'S_N2': 1.4901155163519723e-13, + 'H2O': 0.0} + ''' _N_ins = 0 _N_outs = 0 # _ins_size_is_fixed = True @@ -445,16 +479,7 @@ class BatchExperiment(SanUnit): def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='WasteStream', model=None, isdynamic=True, exogenous_vars=(), **kwargs): - ''' - A batch reactor in experimental settings. - - Parameters - ---------- - model : :class:`CompiledProcesses`, optional - Process model that describes the dynamics of state variables. - The `state` of the batch reactor is entirely determined by the - stoichiometry and rate function in this model. - ''' + SanUnit.__init__(self, ID, None, (), thermo, init_with, isdynamic=isdynamic, exogenous_vars=exogenous_vars) self._model = model @@ -530,233 +555,233 @@ def dy_dt(t, QC_ins, QC, dQC_ins): #TODO: add functions for convenient model calibration #%% NOT READY -class SBR(SanUnit): - ''' - Sequential batch reactors operated in parallel. The number of reactors is - determined by operation cycle and influent flowrate. [1]_ - - Parameters - ---------- - ID : str - ID for the reactors. The default is ''. - ins : :class:`WasteStream` - Influent to the reactor. Expected number of influent is 1. - outs : :class:`WasteStream` - Treated effluent and wasted sludge. - surface_area : float, optional - Surface area of the reactor bottom, in [m^2]. The reactor is assumed - to be cylinder. The default is 1500. - height : float, optional - Height of the reactor, in [m]. The default is 4. - operation_cycle : iterable of float, optional - Operation cycle of the SBR, time for each stage specified in [h]. There - are 7 stages: 1 - fill, 2 - fill, 3 - mix, 4 - mix, 5 - settle, 6 - decant, - 7 - desludge. The first 4 stages are modeled as a biological reactor. - The 5th stage is modeled as a 1D N-layer settler. The last 2 stages are - assumed inactive. The default is (0.5, 1.5, 2.0, 0, 1.0, 0.5, 0.1). - aeration : iterable of float and/or :class:`Process`, optional - Aeration settings for the first 4 stages of the operation cycle. Either - specify a targeted dissolved oxygen concentration in [mg O2/L] or provide - a :class:`Process` object to represent aeration, or None for no aeration. - The default is (None, None, None, 2.0). - DO_ID : str, optional - The :class:`Component` ID for dissolved oxygen, only relevant when the - reactor is aerated. The default is 'S_O2'. - suspended_growth_model : :class:`Processes`, optional - The suspended growth biokinetic model. The default is None. - N_layer : int, optional - The number of layers to model settling. The default is 10. - pumped_flow : float, optional - Designed effluent flowrate, in [m^3/d]. The default is None. - underflow : float, optional - Designed wasted activated sludge flowrate, in [m^3/d]. The default is None. - X_threshold : float, optional - Threshold suspended solid concentration, in [g/m^3]. The default is 3000. - v_max : float, optional - Maximum theoretical (i.e. Vesilind) settling velocity, in [m/d]. The - default is 474. - v_max_practical : float, optional - Maximum practical settling velocity, in [m/d]. The default is 250. - rh : float, optional - Hindered zone settling parameter in the double-exponential settling velocity - function, in [m^3/g]. The default is 5.76e-4. - rp : float, optional - Flocculant zone settling parameter in the double-exponential settling velocity - function, in [m^3/g]. The default is 2.86e-3. - fns : float, optional - Non-settleable fraction of the suspended solids, dimensionless. Must be within - [0, 1]. The default is 2.28e-3. - cache_state : bool, optional - Whether to store volume and composition of retained sludge in the tank from - most recent run. The default is True. - - References - ---------- - .. [1] Takács, I.; Patry, G. G.; Nolasco, D. A Dynamic Model of the Clarification - -Thickening Process. Water Res. 1991, 25 (10), 1263–1271. - https://doi.org/10.1016/0043-1354(91)90066-Y. - - ''' - - _N_ins = 1 - _N_outs = 2 - - def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='WasteStream', - surface_area=1500, height=4, - operation_cycle=(0.5, 1.5, 2.0, 0, 1.0, 0.5, 0.1), - aeration=(None, None, None, 2.0), DO_ID='S_O2', - suspended_growth_model=None, N_layer=10, - pumped_flow=None, underflow=None, - X_threshold=3000, v_max=474, v_max_practical=250, - rh=5.76e-4, rp=2.86e-3, fns=2.28e-3, - cache_state=True, **kwargs): - SanUnit.__init__(self, ID, ins, outs, thermo, init_with) - - self._V = surface_area * height - self._A = surface_area - self._h = height - self._operation_cycle = operation_cycle - self._aeration = aeration - self._DO_ID = DO_ID - self._model = suspended_growth_model - self._N_layer = N_layer - self._Q_e = pumped_flow - self._Q_WAS = underflow - self._X_t = X_threshold - self._v_max = v_max - self._v_max_p = v_max_practical - self._rh = rh - self._rp = rp - self._fns = fns - self._cache_state = cache_state - - for attr, value in kwargs.items(): - setattr(self, attr, value) - self._init_Vas = None - self._init_Cas = None - self._dynamic_composition = None - - - @property - def operation_cycle(self): - return dict(zip(('fill_1', 'fill_2', 'mix_1', 'mix_2', 'settle', 'decant', 'desludge'), - self._operation_cycle)) - @property - def total_cycle_time(self): - return sum(self._operation_cycle) - - @property - def aeration(self): - return dict(zip(('fill_1', 'fill_2', 'mix_1', 'mix_2'), - self._aeration[:4])) - - @property - def C_t(self): - if self._dynamic_composition: - return pd.DataFrame(self._dynamic_composition, - columns = ['Time[d]'] + list(self.components.IDs)) - else: return None - - def _run(self, cache_state=True): - if self._model is None: - raise RuntimeError(f'{self.ID} was initialized without a suspended growth model.') - else: - isa = isinstance - inf = self.ins[0] - Q_in = inf.get_total_flow('m3/d') - eff, sludge = self.outs - eff.copy_like(inf) - sludge.copy_like(inf) - C_in = inf.mass / inf.F_vol * 1e3 # concentrations in g/m3 - cmps = self.components - C = list(symbols(cmps.IDs)) - if self._init_Vas is not None: - V_0 = self._init_Vas - C_0 = self._init_Cas - else: - V_0 = 0 - C_0 = C_in - n = self._N_layer - if self._aeration.count(None) == len(self._aeration): - Vmax = self._V - hj = self._h/n - else: - Vmax = self._V*0.75 - hj = self._h*0.75/n - - # ********fill and mix/aerate stages*********** - T_fill = (Vmax - V_0)/Q_in # maximum total fill time in day - T = [t/24 for t in self._operation_cycle] # operation cycle in day - if T_fill <= T[0]: - schedule = [T_fill, T[0]-T_fill] + T[1:4] - aer = [self._aeration[0], self._aeration[0]] + list(self._aeration[1:4]) - fill = [True] + [False]*4 - V_total = Vmax - elif T_fill <= T[0]+T[1]: - schedule = [T[0], T_fill-T[0], T[0]+T[1]-T_fill] + T[2:4] - aer = list(self._aeration[:2]) + [self._aeration[1]] + list(self._aeration[2:4]) - fill = [True]*2 + [False]*3 - V_total = Vmax - else: - schedule = T[:4] - aer = list(self._aeration[:4]) - fill = [True]*2 + [False]*2 - V_total = Q_in*(T[0]+T[1])+V_0 - hj = V_total/self._V*self._h/n - - for i in range(1, len(schedule)): - if fill[-i] == fill[-i-1] and aer[-i] == aer[-i-1]: - schedule[-i-1] += schedule[-i] - schedule[-i] = 0 - - t_arr = np.array([]) - y_mat = np.ndarray([]) - for i in range(len(schedule)): - if schedule[i] > 0: - dC_dt, J_func = self._compile_dC_dt(V_0, Q_in, C_in, C, fill[i], aer[i]) - if isa(aer[i], (float, int)): C_0[cmps.index(self._DO_ID)] = aer[i] - sol = solve_ivp(dC_dt, (0, schedule[i]), C_0, method='BDF', jac=J_func) - C_0 = sol.y.transpose()[-1] - V_0 += Q_in * schedule[i] * fill[i] - t_arr = np.concatenate((t_arr, sol.t + t_arr[-1])) - y_mat = np.hstack((y_mat, sol.y)) - self._dynamic_composition = np.vstack((t_arr, y_mat)).transpose() - - # *********settle, decant, desludge********** - eff.set_flow(C_0*eff.F_vol, 'g/hr', self.components.IDs) - X_0 = eff.get_TSS() - X_min = X_0 * self._fns - T_settle = T[4] - def dX_dt(t, X): - VX = [_settling_flux(x, self._v_max, self._v_max_p, X_min, self._rh, self._rp) for x in X] - J = [VX[j] if X[j+1] <= self._X_t else min(VX[j], VX[j+1]) for j in range(n-1)] - settle_out = np.array(J + [0]) - settle_in = np.array([0] + J) - dXdt = (settle_in - settle_out)/hj - return dXdt - sol = solve_ivp(dX_dt, (0, T_settle), np.ones(n)*X_0) - X = sol.y.transpose()[-1] - - V_eff = min(T[5]*self._Q_e, V_total*(n-1)/n) - n_eff = V_eff/V_total - w_eff = np.array([1]*floor(n_eff)+[n_eff-floor(n_eff)]) - X_eff = np.average(X[:ceil(n_eff)], weights=w_eff) - eff_mass_flow = (X_eff/X_0*cmps.x + (1-cmps.x))*C_0*V_eff/T[5] - eff.set_flow(eff_mass_flow, 'g/d', cmps.IDs) - - V_was = min(T[6]*self._Q_WAS, V_total-V_eff) - X_as = (V_total*X_0 - V_eff*X_eff) / (V_total-V_eff) - C_as = (X_as/X_0*cmps.x + (1-cmps.x))*C_0 - was_mass_flow = C_as*V_was/T[6] - sludge.set_flow(was_mass_flow, 'g/d', cmps.IDs) - - if self._cache_state: - self._init_Vas = V_total - V_eff - V_was - self._init_Cas = C_as - - - def _design(self): - pass +# class SBR(SanUnit): +# ''' +# Sequential batch reactors operated in parallel. The number of reactors is +# determined by operation cycle and influent flowrate. [1]_ + +# Parameters +# ---------- +# ID : str +# ID for the reactors. The default is ''. +# ins : :class:`WasteStream` +# Influent to the reactor. Expected number of influent is 1. +# outs : :class:`WasteStream` +# Treated effluent and wasted sludge. +# surface_area : float, optional +# Surface area of the reactor bottom, in [m^2]. The reactor is assumed +# to be cylinder. The default is 1500. +# height : float, optional +# Height of the reactor, in [m]. The default is 4. +# operation_cycle : iterable of float, optional +# Operation cycle of the SBR, time for each stage specified in [h]. There +# are 7 stages: 1 - fill, 2 - fill, 3 - mix, 4 - mix, 5 - settle, 6 - decant, +# 7 - desludge. The first 4 stages are modeled as a biological reactor. +# The 5th stage is modeled as a 1D N-layer settler. The last 2 stages are +# assumed inactive. The default is (0.5, 1.5, 2.0, 0, 1.0, 0.5, 0.1). +# aeration : iterable of float and/or :class:`Process`, optional +# Aeration settings for the first 4 stages of the operation cycle. Either +# specify a targeted dissolved oxygen concentration in [mg O2/L] or provide +# a :class:`Process` object to represent aeration, or None for no aeration. +# The default is (None, None, None, 2.0). +# DO_ID : str, optional +# The :class:`Component` ID for dissolved oxygen, only relevant when the +# reactor is aerated. The default is 'S_O2'. +# suspended_growth_model : :class:`Processes`, optional +# The suspended growth biokinetic model. The default is None. +# N_layer : int, optional +# The number of layers to model settling. The default is 10. +# pumped_flow : float, optional +# Designed effluent flowrate, in [m^3/d]. The default is None. +# underflow : float, optional +# Designed wasted activated sludge flowrate, in [m^3/d]. The default is None. +# X_threshold : float, optional +# Threshold suspended solid concentration, in [g/m^3]. The default is 3000. +# v_max : float, optional +# Maximum theoretical (i.e. Vesilind) settling velocity, in [m/d]. The +# default is 474. +# v_max_practical : float, optional +# Maximum practical settling velocity, in [m/d]. The default is 250. +# rh : float, optional +# Hindered zone settling parameter in the double-exponential settling velocity +# function, in [m^3/g]. The default is 5.76e-4. +# rp : float, optional +# Flocculant zone settling parameter in the double-exponential settling velocity +# function, in [m^3/g]. The default is 2.86e-3. +# fns : float, optional +# Non-settleable fraction of the suspended solids, dimensionless. Must be within +# [0, 1]. The default is 2.28e-3. +# cache_state : bool, optional +# Whether to store volume and composition of retained sludge in the tank from +# most recent run. The default is True. + +# References +# ---------- +# .. [1] Takács, I.; Patry, G. G.; Nolasco, D. A Dynamic Model of the Clarification +# -Thickening Process. Water Res. 1991, 25 (10), 1263–1271. +# https://doi.org/10.1016/0043-1354(91)90066-Y. + +# ''' + +# _N_ins = 1 +# _N_outs = 2 + +# def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='WasteStream', +# surface_area=1500, height=4, +# operation_cycle=(0.5, 1.5, 2.0, 0, 1.0, 0.5, 0.1), +# aeration=(None, None, None, 2.0), DO_ID='S_O2', +# suspended_growth_model=None, N_layer=10, +# pumped_flow=None, underflow=None, +# X_threshold=3000, v_max=474, v_max_practical=250, +# rh=5.76e-4, rp=2.86e-3, fns=2.28e-3, +# cache_state=True, **kwargs): +# SanUnit.__init__(self, ID, ins, outs, thermo, init_with) + +# self._V = surface_area * height +# self._A = surface_area +# self._h = height +# self._operation_cycle = operation_cycle +# self._aeration = aeration +# self._DO_ID = DO_ID +# self._model = suspended_growth_model +# self._N_layer = N_layer +# self._Q_e = pumped_flow +# self._Q_WAS = underflow +# self._X_t = X_threshold +# self._v_max = v_max +# self._v_max_p = v_max_practical +# self._rh = rh +# self._rp = rp +# self._fns = fns +# self._cache_state = cache_state + +# for attr, value in kwargs.items(): +# setattr(self, attr, value) +# self._init_Vas = None +# self._init_Cas = None +# self._dynamic_composition = None + + +# @property +# def operation_cycle(self): +# return dict(zip(('fill_1', 'fill_2', 'mix_1', 'mix_2', 'settle', 'decant', 'desludge'), +# self._operation_cycle)) +# @property +# def total_cycle_time(self): +# return sum(self._operation_cycle) + +# @property +# def aeration(self): +# return dict(zip(('fill_1', 'fill_2', 'mix_1', 'mix_2'), +# self._aeration[:4])) + +# @property +# def C_t(self): +# if self._dynamic_composition: +# return pd.DataFrame(self._dynamic_composition, +# columns = ['Time[d]'] + list(self.components.IDs)) +# else: return None + +# def _run(self, cache_state=True): +# if self._model is None: +# raise RuntimeError(f'{self.ID} was initialized without a suspended growth model.') +# else: +# isa = isinstance +# inf = self.ins[0] +# Q_in = inf.get_total_flow('m3/d') +# eff, sludge = self.outs +# eff.copy_like(inf) +# sludge.copy_like(inf) +# C_in = inf.mass / inf.F_vol * 1e3 # concentrations in g/m3 +# cmps = self.components +# C = list(symbols(cmps.IDs)) +# if self._init_Vas is not None: +# V_0 = self._init_Vas +# C_0 = self._init_Cas +# else: +# V_0 = 0 +# C_0 = C_in +# n = self._N_layer +# if self._aeration.count(None) == len(self._aeration): +# Vmax = self._V +# hj = self._h/n +# else: +# Vmax = self._V*0.75 +# hj = self._h*0.75/n + +# # ********fill and mix/aerate stages*********** +# T_fill = (Vmax - V_0)/Q_in # maximum total fill time in day +# T = [t/24 for t in self._operation_cycle] # operation cycle in day +# if T_fill <= T[0]: +# schedule = [T_fill, T[0]-T_fill] + T[1:4] +# aer = [self._aeration[0], self._aeration[0]] + list(self._aeration[1:4]) +# fill = [True] + [False]*4 +# V_total = Vmax +# elif T_fill <= T[0]+T[1]: +# schedule = [T[0], T_fill-T[0], T[0]+T[1]-T_fill] + T[2:4] +# aer = list(self._aeration[:2]) + [self._aeration[1]] + list(self._aeration[2:4]) +# fill = [True]*2 + [False]*3 +# V_total = Vmax +# else: +# schedule = T[:4] +# aer = list(self._aeration[:4]) +# fill = [True]*2 + [False]*2 +# V_total = Q_in*(T[0]+T[1])+V_0 +# hj = V_total/self._V*self._h/n + +# for i in range(1, len(schedule)): +# if fill[-i] == fill[-i-1] and aer[-i] == aer[-i-1]: +# schedule[-i-1] += schedule[-i] +# schedule[-i] = 0 + +# t_arr = np.array([]) +# y_mat = np.ndarray([]) +# for i in range(len(schedule)): +# if schedule[i] > 0: +# dC_dt, J_func = self._compile_dC_dt(V_0, Q_in, C_in, C, fill[i], aer[i]) +# if isa(aer[i], (float, int)): C_0[cmps.index(self._DO_ID)] = aer[i] +# sol = solve_ivp(dC_dt, (0, schedule[i]), C_0, method='BDF', jac=J_func) +# C_0 = sol.y.transpose()[-1] +# V_0 += Q_in * schedule[i] * fill[i] +# t_arr = np.concatenate((t_arr, sol.t + t_arr[-1])) +# y_mat = np.hstack((y_mat, sol.y)) +# self._dynamic_composition = np.vstack((t_arr, y_mat)).transpose() + +# # *********settle, decant, desludge********** +# eff.set_flow(C_0*eff.F_vol, 'g/hr', self.components.IDs) +# X_0 = eff.get_TSS() +# X_min = X_0 * self._fns +# T_settle = T[4] +# def dX_dt(t, X): +# VX = [_settling_flux(x, self._v_max, self._v_max_p, X_min, self._rh, self._rp) for x in X] +# J = [VX[j] if X[j+1] <= self._X_t else min(VX[j], VX[j+1]) for j in range(n-1)] +# settle_out = np.array(J + [0]) +# settle_in = np.array([0] + J) +# dXdt = (settle_in - settle_out)/hj +# return dXdt +# sol = solve_ivp(dX_dt, (0, T_settle), np.ones(n)*X_0) +# X = sol.y.transpose()[-1] + +# V_eff = min(T[5]*self._Q_e, V_total*(n-1)/n) +# n_eff = V_eff/V_total +# w_eff = np.array([1]*floor(n_eff)+[n_eff-floor(n_eff)]) +# X_eff = np.average(X[:ceil(n_eff)], weights=w_eff) +# eff_mass_flow = (X_eff/X_0*cmps.x + (1-cmps.x))*C_0*V_eff/T[5] +# eff.set_flow(eff_mass_flow, 'g/d', cmps.IDs) + +# V_was = min(T[6]*self._Q_WAS, V_total-V_eff) +# X_as = (V_total*X_0 - V_eff*X_eff) / (V_total-V_eff) +# C_as = (X_as/X_0*cmps.x + (1-cmps.x))*C_0 +# was_mass_flow = C_as*V_was/T[6] +# sludge.set_flow(was_mass_flow, 'g/d', cmps.IDs) + +# if self._cache_state: +# self._init_Vas = V_total - V_eff - V_was +# self._init_Cas = C_as + + +# def _design(self): +# pass # def _compile_dC_dt(self, V0, Qin, Cin, C, fill, aer): # isa = isinstance From d23de50b044bfff548e3010ec015501a0869dba9 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Wed, 23 Oct 2024 10:19:29 -0700 Subject: [PATCH 07/12] try different doctest format --- .../sanunits/_suspended_growth_bioreactor.py | 30 +++++++++---------- 1 file changed, 14 insertions(+), 16 deletions(-) diff --git a/qsdsan/sanunits/_suspended_growth_bioreactor.py b/qsdsan/sanunits/_suspended_growth_bioreactor.py index 0196daf8..4b8cc3e8 100644 --- a/qsdsan/sanunits/_suspended_growth_bioreactor.py +++ b/qsdsan/sanunits/_suspended_growth_bioreactor.py @@ -454,22 +454,20 @@ class BatchExperiment(SanUnit): >>> BE = su.BatchExperiment('BE', model=asm1) >>> BE.set_init_conc(S_S=20, X_BH=500, S_O=8, S_ND=3, S_ALK=84) >>> BE.simulate(t_span=(0,10), method='BDF') - >>> BE.state - {'S_I': 0.0, - 'S_S': 0.9259111676784925, - 'X_I': 0.0, - 'X_S': 446.163046741642, - 'X_BH': 25.664835161702214, - 'X_BA': 0.0, - 'X_P': 39.246207146417575, - 'S_O': -1.7987340043125508e-08, - 'S_NO': -4.085240753452882e-20, - 'S_NH': 2.1230597299994245, - 'S_ND': 1.6652546123923254e-06, - 'X_ND': 36.46897947131222, - 'S_ALK': 85.82051689686057, - 'S_N2': 1.4901155163519723e-13, - 'H2O': 0.0} + >>> for k,v in BE.state.items(): + ... if v != 0: + ... print(f'{k}{" "*(7-len(k))}{v:.2f}') + S_S 0.93 + X_S 446.16 + X_BH 25.66 + X_P 39.25 + S_O -0.00 + S_NO -0.00 + S_NH 2.12 + S_ND 0.00 + X_ND 36.47 + S_ALK 85.82 + S_N2 0.00 ''' _N_ins = 0 From 8c2998eba231c4535e291f614d5ef0bae4e02fe9 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Wed, 23 Oct 2024 10:45:20 -0700 Subject: [PATCH 08/12] Update _adm1_p_extension.tsv --- qsdsan/data/process_data/_adm1_p_extension.tsv | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qsdsan/data/process_data/_adm1_p_extension.tsv b/qsdsan/data/process_data/_adm1_p_extension.tsv index 5b91bb98..5ac6f4d8 100644 --- a/qsdsan/data/process_data/_adm1_p_extension.tsv +++ b/qsdsan/data/process_data/_adm1_p_extension.tsv @@ -1,4 +1,4 @@ - S_su S_aa S_fa S_va S_bu S_pro S_ac S_h2 S_ch4 S_IC S_IN S_IP S_I X_ch X_pr X_li X_su X_aa X_fa X_c4 X_pro X_ac X_h2 X_I X_PHA X_PP X_PAO S_K S_Mg X_MeOH X_MeP + S_su S_aa S_fa S_va S_bu S_pro S_ac S_h2 S_ch4 S_IC S_IN S_IP S_I X_ch X_pr X_li X_su X_aa X_fa X_c4 X_pro X_ac X_h2 X_I X_PHA X_PP X_PAO S_K S_Mg hydrolysis_carbs 1 ? ? ? -1 hydrolysis_proteins 1 ? ? ? -1 hydrolysis_lipids 1-f_fa_li f_fa_li ? ? ? -1 From 93a15355a87775b928c90dd73df281065815682b Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Wed, 23 Oct 2024 10:50:28 -0700 Subject: [PATCH 09/12] temporarily disable `mADM1` on main branch --- qsdsan/processes/__init__.py | 6 +- qsdsan/processes/_madm1.py | 787 ----------------------------------- 2 files changed, 3 insertions(+), 790 deletions(-) delete mode 100644 qsdsan/processes/_madm1.py diff --git a/qsdsan/processes/__init__.py b/qsdsan/processes/__init__.py index 5aacbc5a..ba9beecb 100644 --- a/qsdsan/processes/__init__.py +++ b/qsdsan/processes/__init__.py @@ -63,7 +63,7 @@ def __init__(self): from ._asm2d import * from ._adm1 import * from ._adm1_p_extension import * -from ._madm1 import * +# from ._madm1 import * from ._decay import * from ._kinetic_reaction import * from ._pm2 import * @@ -74,7 +74,7 @@ def __init__(self): _asm2d, _adm1, _adm1_p_extension, - _madm1, + # _madm1, _decay, _kinetic_reaction, _pm2 @@ -86,7 +86,7 @@ def __init__(self): *_asm2d.__all__, *_adm1.__all__, *_adm1_p_extension.__all__, - *_madm1.__all__, + # *_madm1.__all__, *_decay.__all__, *_kinetic_reaction.__all__, *_pm2.__all__, diff --git a/qsdsan/processes/_madm1.py b/qsdsan/processes/_madm1.py deleted file mode 100644 index 02ff1b77..00000000 --- a/qsdsan/processes/_madm1.py +++ /dev/null @@ -1,787 +0,0 @@ -# -*- coding: utf-8 -*- -''' -QSDsan: Quantitative Sustainable Design for sanitation and resource recovery systems - -This module is developed by: - 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 chemicals.elements import molecular_weight as get_mw -from qsdsan import Component, Components, Process, Processes, CompiledProcesses -import numpy as np, qsdsan.processes as pc, qsdsan as qs -from qsdsan.utils import ospath, data_path -from qsdsan.processes._adm1 import ( - R, - create_adm1_cmps, - ADM1, - mass2mol_conversion, - T_correction_factor, - substr_inhibit, - non_compet_inhibit, - Hill_inhibit - ) -# from scipy.optimize import brenth -# from warnings import warn - - -__all__ = ('create_madm1_cmps', 'ModifiedADM1') - -_path = ospath.join(data_path, 'process_data/_madm1.tsv') - -#%% components -# C_mw = get_mw({'C':1}) -N_mw = get_mw({'N':1}) -P_mw = get_mw({'P':1}) -S_mw = get_mw({'S':1}) -Fe_mw = get_mw({'Fe':1}) -O_mw = get_mw({'O':1}) - -def create_madm1_cmps(set_thermo=True, ASF_L=0.31, ASF_H=1.2): - ''' - Create a set of components for the modified ADM1. - - Parameters - ---------- - set_thermo : bool, optional - Whether to set thermo with the returned set of components. The default is True. - ASF_L : float, optional - Active site factor for X_HFO_L [mol P sites/mol Fe]. The default is 0.31. - ASF_H : float, optional - Active site factor for X_HFO_H [mol P sites/mol Fe]. The default is 1.2. - - Returns - ------- - cmps_madm1 : class:`CompiledComponents` - - ''' - - # Components from the original ADM1 - # ********************************* - _cmps = create_adm1_cmps(False) - S_aa = _cmps.S_aa - X_pr = _cmps.X_pr - S_aa.i_C = X_pr.i_C = 0.36890 - S_aa.i_N = X_pr.i_N = 0.11065 - S_aa.i_P = X_pr.i_P = 0. - S_aa.i_mass = X_pr.i_mass = 1/1.35566 - - S_fa = _cmps.S_fa - S_fa.formula = 'C25H52O3' - S_bu = _cmps.S_bu - S_bu.formula = 'C4H8O2' - S_pro = _cmps.S_pro - S_pro.formula = 'C3H6O2' - S_ac = _cmps.S_ac - S_ac.formula = 'C2H4O2' - - S_I = _cmps.S_I - X_I = _cmps.X_I - S_I.i_C = X_I.i_C = 0.36178 - S_I.i_N = X_I.i_N = 0.06003 - S_I.i_P = X_I.i_P = 0.00649 - S_I.i_mass = X_I.i_mass = 1/1.54100 - - X_ch = _cmps.X_ch - X_ch.formula = 'C24H48O24' - # _cmps.X_li.formula = 'C64H119O7.5P' - X_li = X_pr.copy('X_li') - X_li.i_C = 0.26311 - X_li.i_N = 0. - X_li.i_P = 0.01067 - X_li.i_mass = 1/2.81254 - - adm1_biomass = (_cmps.X_su, _cmps.X_aa, _cmps.X_fa, _cmps.X_c4, _cmps.X_pro, _cmps.X_ac, _cmps.X_h2) - for bio in adm1_biomass: - # bio.formula = 'C5H7O2NP0.113' - bio.i_C = 0.36612 - bio.i_N = 0.08615 - bio.i_P = 0.02154 - bio.i_mass = 1/1.39300 - - # P related components from ASM2d - # ******************************* - asm_cmps = pc.create_asm2d_cmps(False) - X_PHA = asm_cmps.X_PHA - X_PHA.formula = '(C2H4O)n' - # X_PHA.i_C = 0.3 - # X_PHA.i_mass = 0.55 - - X_PAO = _cmps.X_su.copy('X_PAO') - X_PAO.description = 'Phosphorus-accumulating organism biomass' - - # Additional components for P, S, Fe extensions - # ********************************************* - S_IP = asm_cmps.S_PO4.copy('S_IP') - - ion_properties = dict( - particle_size='Soluble', - degradability='Undegradable', - organic=False) - S_K = Component.from_chemical('S_K', chemical='K+', description='Potassium ion', - measured_as='K', **ion_properties) - S_Mg = Component.from_chemical('S_Mg', chemical='Mg2+', description='Magnesium ion', - measured_as='Mg',**ion_properties) - S_SO4 = Component.from_chemical('S_SO4', chemical='SO4-2', description='Sulfate', - measured_as='S', **ion_properties) - S_IS = Component.from_chemical('S_IS', chemical='H2S', - description='Hydrogen sulfide', - measured_as='COD', - particle_size='Soluble', - degradability='Undegradable', - organic=False) - - X_hSRB = X_PAO.copy('X_hSRB') - X_hSRB.description = 'sulfate-reducing biomass, utilizing H2' - X_aSRB = X_PAO.copy('X_aSRB') - X_aSRB.description = 'sulfate-reducing biomass, utilizing acetate' - X_pSRB = X_PAO.copy('X_pSRB') - X_pSRB.description = 'sulfate-reducing biomass, utilizing propionate' - X_c4SRB = X_PAO.copy('X_c4SRB') - X_c4SRB.description = 'sulfate-reducing biomass, utilizing butyrate and valerate' - - S_S0 = Component.from_chemical('S_S0', chemical='S', - description='Elemental sulfur', - measured_as='COD', - particle_size='Soluble', - degradability='Undegradable', - organic=False) - S_Fe3 = Component.from_chemical('S_Fe3', chemical='Fe3+', description='Iron (III)', - measured_as='Fe',**ion_properties) - S_Fe2 = Component.from_chemical('S_Fe2', chemical='Fe2+', description='Iron (II)', - measured_as='Fe',**ion_properties) - S_Fe2.i_COD = 0.5*O_mw/Fe_mw - S_Fe2.measured_as = 'COD' - - # Multiple mineral precipitation - # ****************************** - mineral_properties = dict( - particle_size='Particulate', - degradability='Undegradable', - organic=False) - - X_HFO_H = Component('X_HFO_H', formula='FeO(OH)', - description='Hydrous ferric oxide with high number of active sites', - measured_as='Fe',**mineral_properties) - X_HFO_L = X_HFO_H.copy('X_HFO_L') - X_HFO_L.description = 'Hydrous ferric oxide with low number of active sites' - - X_HFO_old = X_HFO_H.copy('X_HFO_old') - X_HFO_old.description = 'Inactive hydrous ferric oxide' - - X_HFO_HP = Component('X_HFO_HP', formula=f'FeO(OH)P{ASF_H}', - description='X_HFO_H with phosphorus-bounded adsorption sites', - measured_as='Fe', **mineral_properties) - X_HFO_HP_old = X_HFO_HP.copy('X_HFO_HP_old') - X_HFO_HP_old.description = 'Old ' + X_HFO_HP.description - - X_HFO_LP = Component('X_HFO_LP', formula=f'FeO(OH)P{ASF_L}', - description='X_HFO_L with phosphorus-bounded adsorption sites', - measured_as='Fe', **mineral_properties) - X_HFO_LP_old = X_HFO_LP.copy('X_HFO_LP_old') - X_HFO_LP_old.description = 'Old ' + X_HFO_LP.description - - X_CCM = Component.from_chemical('X_CCM', chemical='calcite', description='Calcite', **mineral_properties) - X_ACC = Component.from_chemical('X_ACC', chemical='aragonite', description='Aragonite', **mineral_properties) - X_ACP = Component.from_chemical('X_ACP', chemical='Ca3(PO4)2', description='Amorphous calcium phosphate', **mineral_properties) - X_HAP = Component.from_chemical('X_HAP', chemical='hydroxylapatite', description='Hydroxylapatite', **mineral_properties) - X_DCPD = Component.from_chemical('X_DCPD', chemical='CaHPO4', description='Dicalcium phosphate', **mineral_properties) - X_OCP = Component('X_OCP', formula='Ca4HP3O12', description='Octacalcium phosphate', **mineral_properties) - X_struv = Component.from_chemical('X_struv', chemical='MgNH4PO4', description='Struvite', **mineral_properties) - X_newb = Component.from_chemical('X_newb', chemical='MgHPO4', description='Newberyite', **mineral_properties) - X_magn = Component.from_chemical('X_magn', chemical='MgCO3', description='Magnesite', **mineral_properties) - X_kstruv = Component('X_kstruv', formula='MgKPO4', description='K-struvite', **mineral_properties) - X_FeS = Component.from_chemical('X_FeS', chemical='FeS', description='Iron sulfide', **mineral_properties) - X_Fe3PO42 = Component('X_Fe3PO42', formula='Fe3(PO4)2', description='Ferrous phosphate', **mineral_properties) - X_AlPO4 = Component.from_chemical('X_AlPO4', chemical='AlPO4', description='Aluminum phosphate', **mineral_properties) - - S_Ca = Component.from_chemical('S_Ca', chemical='Ca2+', description='Calsium ion', - measured_as='Ca', **ion_properties) - S_Al = Component.from_chemical('S_Al', chemical='Al3+', description='Aluminum ion', - measured_as='Al', **ion_properties) - S_Na = Component.from_chemical('S_Na', chemical='Na+', description='Sodium ion', - measured_as='Na', **ion_properties) - S_Cl = Component.from_chemical('S_Cl', chemical='Cl-', description='Chloride', - measured_as='Cl', **ion_properties) - - cmps_madm1 = Components([_cmps.S_su, S_aa, S_fa, _cmps.S_va, S_bu, - S_pro, S_ac, _cmps.S_h2, _cmps.S_ch4, - _cmps.S_IC, _cmps.S_IN, S_IP, S_I, - X_ch, X_pr, X_li, *adm1_biomass, X_I, - X_PHA, asm_cmps.X_PP, X_PAO, S_K, S_Mg, - S_SO4, S_IS, X_hSRB, X_aSRB, X_pSRB, X_c4SRB, - S_S0, S_Fe3, S_Fe2, X_HFO_H, X_HFO_L, X_HFO_old, - X_HFO_HP, X_HFO_LP, X_HFO_HP_old, X_HFO_LP_old, - S_Ca, S_Al, X_CCM, X_ACC, X_ACP, X_HAP, X_DCPD, - X_OCP, X_struv, X_newb, X_magn, X_kstruv, X_FeS, - X_Fe3PO42, X_AlPO4, - S_Na, S_Cl, _cmps.H2O]) - cmps_madm1.default_compile() - - if set_thermo: qs.set_thermo(cmps_madm1) - return cmps_madm1 - -#%% rate functions - -# https://wiki.dynamita.com/en/biokinetic_process_models#chemical-phosphorus-removal-with-metal-salts-addition-iron-or-aluminium - -# ============================================================================= -# state_variable_indices = { -# 'S_su': 0, 'S_aa': 1, 'S_fa': 2, 'S_va': 3, 'S_bu': 4, 'S_pro': 5, 'S_ac': 6, 'S_h2': 7, -# 'S_ch4': 8, 'S_IC': 9, 'S_IN': 10, 'S_IP': 11, 'S_I': 12, -# 'X_ch': 13, 'X_pr': 14, 'X_li': 15, -# 'X_su': 16, 'X_aa': 17, 'X_fa': 18, 'X_c4': 19, 'X_pro': 20, 'X_ac': 21, 'X_h2': 22, 'X_I': 23, -# 'X_PHA': 24, 'X_PP': 25, 'X_PAO': 26, 'S_K': 27, 'S_Mg': 28, -# 'S_SO4': 29, 'S_IS': 30, 'X_hSRB': 31, 'X_aSRB': 32, 'X_pSRB': 33, 'X_c4SRB': 34, -# 'S_S0': 35, 'S_Fe3': 36, 'S_Fe2': 37, -# 'X_HFO_H': 38, 'X_HFO_L': 39, 'X_HFO_old': 40, 'X_HFO_HP': 41, 'X_HFO_LP': 42, 'X_HFO_HP_old': 43, 'X_HFO_LP_old': 44, -# 'S_Ca': 45, 'S_Al': 46, -# 'X_CCM': 47, 'X_ACC': 48, 'X_ACP': 49, 'X_HAP': 50, 'X_DCPD': 51, 'X_OCP': 52, -# 'X_struv': 53, 'X_newb': 54, 'X_magn': 55, 'X_kstruv': 56, -# 'X_FeS': 57, 'X_Fe3PO42': 58, -# 'X_AlPO4': 59, -# 'S_Na': 60, 'S_Cl': 61, 'H2O': 62 -# } -# ============================================================================= - -def calc_pH(): - pass - -def calc_biogas(): - pass - -def pcm(): - pass - -def saturation_index(): - pass - -rhos = np.zeros(38+8+13+4) # 38 biological + 8 chemical P removal by HFO + 13 MMP + 4 gas transfer -Cs = np.empty(38+8) -sum_stoichios = np.array([2, 2, 5, 9, 3, 8, 3, 3, 2, 3, 2, 2]) - -def rhos_madm1(state_arr, params, T_op): - ks = params['rate_constants'] - Ks = params['half_sat_coeffs'] - K_PP = params['K_PP'] - K_so4 = params['K_so4'] - cmps = params['components'] - # n = len(cmps) - pH_LLs, pH_ULs = params['pH_limits'] - KS_IN = params['KS_IN'] - KS_IP = params['KS_IP'] - KI_nh3 = params['KI_nh3'] - KIs_h2 = params['KIs_h2'] - KIs_h2s = params['KIs_h2s'] - KHb = params['K_H_base'] - Kab = params['Ka_base'] - KH_dH = params['K_H_dH'] - Ka_dH = params['Ka_dH'] - kLa = params['kLa'] - k_cryst = params['k_cryst'] - n_cryst = params['n_cryst'] - Kspb = params['Ksp_base'] - Ksp_dH = params['Ksp_dH'] - T_base = params['T_base'] - - Cs[:7] = state_arr[13:20] # original ADM1 processes - Cs[7:11] = state_arr[19:23] - Cs[11:18] = state_arr[16:23] - Cs[18:23] = X_PAO = state_arr[26] # P extension processes - Cs[23:25] = X_PP, X_PHA = state_arr[[25,24]] - Cs[25:27] = state_arr[31] # S extension processes - Cs[27:29] = state_arr[32] - Cs[29:31] = state_arr[33] - Cs[31:34] = state_arr[34] - Cs[34:36] = Cs[36:38] = Cs[38:40] = Cs[40:42] = state_arr[38:40] # Fe extension processes + HFO module - Cs[42:44] = Cs[44:46] = state_arr[41:43] - - rhos[:46] = ks * Cs - primary_substrates = state_arr[:8] - - rhos[3:11] *= substr_inhibit(primary_substrates, Ks[:8]) - c4 = primary_substrates[[3,4]] - if sum(c4) > 0: rhos[[6,7]] *= c4/sum(c4) - - vfas = primary_substrates[3:7] - rhos[18:22] *= substr_inhibit(vfas, Ks[8]) - if sum(vfas) > 0: rhos[18:22] *= vfas/sum(vfas) - if X_PAO > 0: rhos[18:22] *= substr_inhibit(X_PP/X_PAO, K_PP) - - srb_subs = np.flip(primary_substrates[3:]) - S_SO4, S_IS = state_arr[29:31] - rhos[[25,27,29,31,32]] *= substr_inhibit(srb_subs, Ks[9:13]) * substr_inhibit(S_SO4, K_so4) - if sum(srb_subs[-2:]) > 0: rhos[[31,32]] *= srb_subs[-2:]/sum(srb_subs[-2:]) - - #!!! why divide by 16 or 64? - S_h2 = primary_substrates[-1] - rhos[34:36] *= S_h2 / 16 - rhos[36:38] *= S_IS / 64 - - KPbind, KPdiss = Ks[-2:] - S_IP = state_arr[11] - rhos[40:42] *= substr_inhibit(S_IP, KPbind) - rhos[44:46] *= non_compet_inhibit(S_IP, KPdiss) - - # inhibition factors - # ****************** - unit_conversion = mass2mol_conversion(cmps) - if T_op == T_base: - Ka = Kab - KH = KHb / unit_conversion[[7,8,9,30]] - Ksp = Kspb - else: - T_temp = params.pop('T_op', None) - if T_op == T_temp: - params['T_op'] = T_op - Ka = params['Ka'] - KH = params['KH'] - Ksp = params['Ksp'] - else: - params['T_op'] = T_op - Ka = params['Ka'] = Kab * T_correction_factor(T_base, T_op, Ka_dH) - KH = params['KH'] = KHb * T_correction_factor(T_base, T_op, KH_dH) / unit_conversion[[7,8,9,30]] - Ksp = params['Ksp'] = Kspb * T_correction_factor(T_base, T_op, Ksp_dH) - - S_IN, S_IP = state_arr[[10,11]] - I_nutrients = substr_inhibit(S_IN, KS_IN) * substr_inhibit(S_IP, KS_IP) - rhos[3:11] *= I_nutrients - rhos[[25,27,29,31,32]] *= I_nutrients - -# ============================================================================= -# !!! place holder for PCM (speciation) -# ============================================================================= - pH, nh3, co2, acts = pcm(state_arr, params) - Is_pH = Hill_inhibit(10**(-pH), pH_ULs, pH_LLs) - rhos[3:9] *= Is_pH[0] - rhos[9:11] *= Is_pH[1:3] - rhos[[25,27]] *= Is_pH[3:5] - rhos[[29,31,32]] *= Is_pH[-1] - - Is_h2 = non_compet_inhibit(S_h2, KIs_h2) - rhos[5:9] *= Is_h2 - Inh3 = non_compet_inhibit(nh3, KI_nh3) - rhos[9] *= Inh3 - - Z_h2s = calc_biogas() # should be a function of pH, like co2 and nh3 - Is_h2s = non_compet_inhibit(Z_h2s, KIs_h2s) - rhos[6:11] *= Is_h2s[:5] - rhos[[25,27,29,31,32]] *= Is_h2s[5:] - - # multiple mineral precipitation - # ****************************** - SIs = np.maximum(1.0, saturation_index(acts, Ksp)) # should be an array - rhos[46:59] = k_cryst * state_arr[47:60] * (SIs**(1/sum_stoichios) - 1)**n_cryst - - # gas transfer - # ************ - biogas_S = state_arr[[7,8,9,30]].copy() - biogas_S[2] = co2 / unit_conversion[9] - biogas_S[3] = Z_h2s / unit_conversion[30] - biogas_p = R * T_op * state_arr[63:67] - rhos[-4:] = kLa * (biogas_S - KH * biogas_p) - - return rhos - -#%% modified ADM1 class -_load_components = settings.get_default_chemicals - -def fun(q_aging_H=450.0, q_aging_L=0.1, q_Pcoprec=360, q_Pbinding=0.3, q_diss_H=36.0, q_diss_L=36.0, - K_Pbind=37.2, K_Pdiss=0.93): - ''' - - - Parameters - ---------- - - Returns - ------- - None. - - ''' - pass - -@chemicals_user -class ModifiedADM1(CompiledProcesses): - """ - Modified Anaerobic Digestion Model no.1 [1]_, [2]_, [3]_ - - Parameters - ---------- - f_ch_xb : float, optional - Fraction of carbohydrates as biomass decay product. The default is 0.275. - f_pr_xb : flaot, optional - Fraction of proteins as biomass decay product. The default is 0.275. - f_li_xb : float, optional - Fraction of lipids as biomass decay product. The default is 0.35. - f_xI_xb : float, optional - Fraction of inert particulates as biomass decay product. The default is 0.1. - f_va_pha : float, optional - Fraction of valerate as PHA lysis product. The default is 0.1. - f_bu_pha : float, optional - Fraction of butyrate as PHA lysis product. The default is 0.1. - f_pro_pha : float, optional - Fraction of propionate as PHA lysis product. The default is 0.4. - Y_PO4 : float, optional - Poly-phosphorus (PP) required for PHA storage [kg P/kg COD]. The default is 0.4. - Y_hSRB : float, optional - Sulfide-reducing biomass (SRB) yield of hydrogen uptake [kg COD/kg COD]. - The default is 0.05. - Y_aSRB : float, optional - SRB yield of acetate uptake [kg COD/kg COD]. The default is 0.05. - Y_pSRB : float, optional - SRB yield of propionate uptake [kg COD/kg COD]. The default is 0.04. - Y_c4SRB : float, optional - SRB yield of butyrate or valerate uptake [kg COD/kg COD]. - The default is 0.06. - q_pha : float, optional - Maximum specific rate constant for PHA storage by phosphorus-accumulating - organisms (PAOs) [d^(-1)]. The default is 3.0. - b_pao : float, optional - PAO lysis rate constant [d^(-1)]. The default is 0.2. - b_pp : float, optional - PP lysis rate constant [d^(-1)]. The default is 0.2. - b_pha : float, optional - PHA lysis rate constant [d^(-1)]. The default is 0.2. - K_A : float, optional - Substrate half saturation coefficient for PHA storage [kg COD/m3]. - The default is 4e-3. - K_PP : float, optional - PP half saturation coefficient for PHA storage [kg P (X_PP)/kg COD (X_PHA)]. - The default is 0.01. - k_hSRB : float, optional - Maximum specific growth rate constant of hydrogen-uptaking SRB [d^(-1)]. - The default is 41.125. - k_aSRB : float, optional - Maximum specific growth rate constant of acetate-uptaking SRB [d^(-1)]. - The default is 10.. - k_pSRB : float, optional - Maximum specific growth rate constant of propionate-uptaking SRB [d^(-1)]. - The default is 16.25. - k_c4SRB : float, optional - Maximum specific growth rate constant of butyrate- or valerate-uptaking - SRB [d^(-1)]. The default is 23. - b_hSRB : float, optional - Hydrogen-uptaking SRB decay rate constant [d^(-1)]. The default is 0.02. - b_aSRB : float, optional - Acetate-uptaking SRB decay rate constant [d^(-1)]. The default is 0.02. - b_pSRB : float, optional - Propionate-uptaking SRB decay rate constant [d^(-1)]. The default is 0.02. - b_c4SRB : float, optional - Butyrate- or valerate-uptaking SRB decay rate constant [d^(-1)]. - The default is 0.02. - K_hSRB : float, optional - Substrate half saturation coefficient of hydrogen uptake by SRB - [kg COD/m3]. The default is 5.96e-6. - K_aSRB : float, optional - Substrate half saturation coefficient of acetate uptake by SRB - [kg COD/m3]. The default is 0.176. - K_pSRB : float, optional - Substrate half saturation coefficient of propionate uptake by SRB - [kg COD/m3]. The default is 0.088. - K_c4SRB : float, optional - Substrate half saturation coefficient of butyrate or valerate uptake by - SRB [kg COD/m3]. The default is 0.1739. - K_so4_hSRB : float, optional - Sulfate half saturation coefficient of SRB uptaking hydrogen [kg S/m3]. - The default is 3.335e-3. - K_so4_aSRB : float, optional - Sulfate half saturation coefficient of SRB uptaking acetate [kg S/m3]. - The default is 6.413e-3. - K_so4_pSRB : float, optional - Sulfate half saturation coefficient of SRB uptaking propionate [kg S/m3]. - The default is 6.413e-3. - K_so4_c4SRB : float, optional - Sulfate half saturation coefficient of SRB uptaking butyrate or valerate - [kg S/m3]. The default is 6.413e-3. - k_Fe3t2_h2 : float, optional - Fe(3+) reduction rate constant [m3∙kg^(-1) Fe(III)∙d^(-1)] using hydrogen - as electron donor. The default is 1.79e7. - k_Fe3t2_is : float, optional - Fe(3+) reduction rate constant [m3∙kg^(-1) Fe(III)∙d^(-1)] using sulfide - as electron donor. The default is 1.79e7. - KS_IP : float, optional - Inorganic phosphorus (nutrient) inhibition coefficient for soluble - substrate uptake [M]. The default is 2e-5. - q_aging_H : float, optional - Aging rate constant of X_HFO_H and X_HFO_HP [d^(-1)]. The default is 450.0. - q_aging_L : float, optional - Aging rate constant of X_HFO_L and X_HFO_LP [d^(-1)]. The default is 0.1. - q_Pcoprec : float, optional - Rate constant of P binding and coprecipitation on X_HFO_H [d^(-1)]. - The default is 360. - q_Pbinding : float, optional - Rate constant of P binding on X_HFO_L [d^(-1)]. The default is 0.3. - q_diss_H : float, optional - Dissolution rate constant of X_HFO_HP [d^(-1)]. The default is 36.0. - q_diss_L : float, optional - Dissolution rate constant of X_HFO_HP [d^(-1)]. The default is 36.0. - K_Pbind : float, optional - S_IP half saturation coefficient for binding with X_HFO_H or X_HFO_L - [kg P/m3]. The default is 37.2, i.e., 1.20 kmol P/m3. - K_Pdiss : float, optional - S_IP half inhibition coefficient for dissolution of X_HFO_HP or X_HFO_LP - [kg P/m3]. The default is 0.93, i.e., 0.03 kmol P/m3. - KI_h2s_c4 : float, optional - H2S half inhibition coefficient for butyrate or valerate uptake - [kg COD/m3]. The default is 0.481. - KI_h2s_pro : float, optional - H2S half inhibition coefficient for propionate uptake [kg COD/m3]. - The default is 0.481. - KI_h2s_ac : float, optional - H2S half inhibition coefficient for acetate uptake [kg COD/m3]. - The default is 0.460. - KI_h2s_h2 : float, optional - H2S half inhibition coefficient for hydrogen uptake [kg COD/m3]. - The default is 0.400. - KI_h2s_c4SRB : float, optional - H2S half inhibition coefficient for butyrate or valerate uptake by SRB - [kg COD/m3]. The default is 0.520. - KI_h2s_pSRB : float, optional - H2S half inhibition coefficient for propionate uptake by SRB [kg COD/m3]. - The default is 0.520. - KI_h2s_aSRB : float, optional - H2S half inhibition coefficient for acetate uptake by SRB [kg COD/m3]. - The default is 0.499. - KI_h2s_hSRB : float, optional - H2S half inhibition coefficient for hydrogen uptake by SRB [kg COD/m3]. - The default is 0.499. - pH_limits_aa_SRB : 2-tuple, optional - Lower and upper limits of pH inhibition for acetogenosis by SRB, - unitless. The default is (6,7). - pH_limits_ac_SRB : 2-tuple, optional - Lower and upper limits of pH inhibition for acetate uptake by SRB, - unitless. The default is (6,7). - pH_limits_h2_SRB : 2-tuple, optional - Lower and upper limits of pH inhibition for hydrogen uptake by SRB, - unitless. The default is (5,6). - k_cryst : iterable[float], optional - Mineral precipitation rate constants [h^(-1)], following the order of - `ModifiedADM1._precipitates`. The default is - [0.35, 1e-3, 3.0, 1e-3, 2.0, 0.76, 5.0, 1e-3, 1e-3, 1e-3, 1e2, 1e-3, 1e-3]. - n_cryst : iterable[int], optional - The effect orders of mineral precipitation reactions [unitless], following - the order of `ModifiedADM1._precipitates`. The default is - [2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2]. - - - Examples - -------- - ... - - References - ---------- - .. [1] Flores-Alsina, X., Solon, K., Kazadi Mbamba, C., Tait, S., - Gernaey, K. V., Jeppsson, U., Batstone, D. J. (2016). - Modelling phosphorus (P), sulfur (S) and iron (Fe) interactions - for dynamic simulations of anaerobic digestion processes. - Water Research, 95, 370–382. https://doi.org/10.1016/J.WATRES.2016.03.012 - .. [2] Solon, K., Flores-Alsina, X., Kazadi Mbamba, C., Ikumi, D., - Volcke, E. I. P., Vaneeckhaute, C., Ekama, G., Vanrolleghem, P. A., - Batstone, D. J., Gernaey, K. v., Jeppsson, U. (2017). Plant-wide - modelling of phosphorus transformations in wastewater treatment systems: - Impacts of control and operational strategies. Water Research, - 113, 97–110. https://doi.org/10.1016/J.WATRES.2017.02.007 - .. [3] Hauduc, H., Takács, I., Smith, S., Szabo, A., Murthy, S., Daigger, G. T., - Spérandio, M. (2015). A dynamic physicochemical model for chemical phosphorus - removal. Water Research, 73, 157–170. https://doi.org/10.1016/J.WATRES.2014.12.053 - - See Also - -------- - `qsdsan.processes.ADM1 `_ - - """ - - _cmp_dependent_stoichio = ('K_XPP', 'Mg_XPP', - 'MW_S0', 'MW_IS', - 'i_mass_S0', 'i_mass_IS', 'i_mass_Fe2') - _stoichio_params = (*ADM1._stoichio_params[5:], - 'f_ch_xb', 'f_pr_xb', 'f_li_xb', 'f_xI_xb', 'f_sI_xb', - 'f_va_pha', 'f_bu_pha', 'f_pro_pha', 'f_ac_pha', - 'f_is_pro', 'f_is_bu', 'f_is_va', - 'Y_PO4', 'Y_hSRB', 'Y_aSRB', 'Y_pSRB', 'Y_c4SRB', - *_cmp_dependent_stoichio - ) - _kinetic_params = ('rate_constants', 'half_sat_coeffs', 'K_PP', 'K_so4', - 'pH_limits', 'KS_IN', 'KS_IP', 'KI_nh3', 'KIs_h2', 'KIs_h2s' - 'Ka_base', 'Ka_dH', 'K_H_base', 'K_H_dH', 'kLa', - 'k_cryst', 'n_cryst', 'Ksp_base', 'Ksp_dH', - 'T_base', 'components', - # 'root' - ) - _acid_base_pairs = ADM1._acid_base_pairs - _biogas_IDs = (*ADM1._biogas_IDs, 'S_IS') - _biomass_IDs = (*ADM1._biomass_IDs, 'X_PAO', 'X_hSRB', 'X_aSRB', 'X_pSRB', 'X_c4SRB') - _precipitates = ('X_CCM', 'X_ACC', 'X_ACP', 'X_HAP', 'X_DCPD', 'X_OCP', - 'X_struv', 'X_newb', 'X_magn', 'X_kstruv', - 'X_FeS', 'X_Fe3PO42', 'X_AlPO4') - _T_base = 298.15 - _K_H_base = [7.8e-4, 1.4e-3, 3.5e-2, 0.105] # biogas species Henry's Law constant [M/bar] - _K_H_dH = [-4180, -14240, -19410, -19180] # Heat of reaction of liquid-gas transfer of biogas species [J/mol] - - _pKsp_base = [8.48, 8.3, 28.92, 44.333, 18.995, 47.08, - 13.6, 18.175, 7.46, 11.5508, - 2.95, 37.76, 18.2] - _Ksp_dH = [8000, -12000, 54000, 0, 31000, 0, - -22600, -22600, -20000, -22600, - -11000, 5060, 0] - - def __new__(cls, components=None, path=None, - f_ch_xb=0.275, f_pr_xb=0.275, f_li_xb=0.35, f_xI_xb=0.1, - f_fa_li=0.95, f_bu_su=0.1328, f_pro_su=0.2691, f_ac_su=0.4076, - f_va_aa=0.23, f_bu_aa=0.26, f_pro_aa=0.05, f_ac_aa=0.4, - f_ac_fa=0.7, f_pro_va=0.54, f_ac_va=0.31, f_ac_bu=0.8, f_ac_pro=0.57, - Y_su=0.1, Y_aa=0.08, Y_fa=0.06, Y_c4=0.06, Y_pro=0.04, Y_ac=0.05, Y_h2=0.06, - f_va_pha=0.1, f_bu_pha=0.1, f_pro_pha=0.4, - Y_PO4=0.4, Y_hSRB=0.05, Y_aSRB=0.05, Y_pSRB=0.04, Y_c4SRB=0.06, - q_ch_hyd=10, q_pr_hyd=10, q_li_hyd=10, - k_su=30, k_aa=50, k_fa=6, k_c4=20, k_pro=13, k_ac=8, k_h2=35, - K_su=0.5, K_aa=0.3, K_fa=0.4, K_c4=0.2, K_pro=0.1, K_ac=0.15, K_h2=7e-6, - b_su=0.02, b_aa=0.02, b_fa=0.02, b_c4=0.02, b_pro=0.02, b_ac=0.02, b_h2=0.02, - q_pha=3.0, b_pao=0.2, b_pp=0.2, b_pha=0.2, K_A=4e-3, K_PP=0.01, - k_hSRB=41.125, k_aSRB=10., k_pSRB=16.25, k_c4SRB=23, - b_hSRB=0.02, b_aSRB=0.02, b_pSRB=0.02, b_c4SRB=0.02, - K_hSRB=5.96e-6, K_aSRB=0.176, K_pSRB=0.088, K_c4SRB=0.1739, - K_so4_hSRB=1.04e-4*S_mw, K_so4_aSRB=2e-4*S_mw, K_so4_pSRB=2e-4*S_mw, K_so4_c4SRB=2e-4*S_mw, - k_Fe3t2_h2=1e9/Fe_mw, k_Fe3t2_is=1e9/Fe_mw, - q_aging_H=450.0, q_aging_L=0.1, q_Pcoprec=360, q_Pbinding=0.3, q_diss_H=36.0, q_diss_L=36.0, - K_Pbind=37.2, K_Pdiss=0.93, # 1.20 and 0.03 in MATLAB, assuming in kmol-P/m3 ? - KI_h2_fa=5e-6, KI_h2_c4=1e-5, KI_h2_pro=3.5e-6, KI_nh3=1.8e-3, KS_IN=1e-4, KS_IP=2e-5, - KI_h2s_c4=0.481, KI_h2s_pro=0.481, KI_h2s_ac=0.460, KI_h2s_h2=0.400, - KI_h2s_c4SRB=0.520, KI_h2s_pSRB=0.520, KI_h2s_aSRB=0.499, KI_h2s_hSRB=0.499, - pH_limits_aa=(4,5.5), pH_limits_ac=(6,7), pH_limits_h2=(5,6), - pH_limits_aa_SRB=(6,7), pH_limits_ac_SRB=(6,7), pH_limits_h2_SRB=(5,6), - kLa=200, pKa_base=[14, 9.25, 6.35, 4.76, 4.88, 4.82, 4.86], - Ka_dH=[55900, 51965, 7646, 0, 0, 0, 0], - k_cryst=[0.35, 1e-3, 3.0, 1e-3, 2.0, 0.76, 5.0, 1e-3, 1e-3, 1e-3, 1e2, 1e-3, 1e-3], - n_cryst=[2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2], - **kwargs): - - cmps = _load_components(components) - - if not path: path = _path - self = Processes.load_from_file(path, - components=cmps, - conserved_for=('C', 'N', 'P'), - parameters=cls._stoichio_params, - compile=False) - - for i in ('fast_P_binding', 'slow_P_sorption', 'dissolution_HFO_HP', 'dissolution_HFO_LP'): - p = getattr(self, i) - p.ref_component = 'S_IP' - - precipitation = [] - for i in cls._precipitates[:-3]: - new_p = Process('precipitation_%s' % i.lstrip('X_'), - reaction='[?]S_IC + [?]S_IN + [?]S_IP + [?]S_K + [?]S_Mg + [?]S_Ca -> %s' % i, - ref_component=i, - conserved_for=('C', 'N', 'P', 'K', 'Mg', 'Ca'), - parameters=()) - precipitation.append(new_p) - - i_mass_IS = cmps.S_IS.i_mass - i_mass_Fe2 = cmps.S_Fe2.i_mass - FeS_mw = cmps.X_FeS.chem_MW - new_p = Process('precipitation_FeS', - reaction={'S_Fe2': -Fe_mw/FeS_mw/i_mass_Fe2, - 'S_IS': -S_mw/FeS_mw/i_mass_IS, - 'X_FeS': 1}, - ref_component='X_FeS', - conserved_for=()) - precipitation.append(new_p) - - Fe3PO42_mw = cmps.X_Fe3PO42.chem_MW - new_p = Process('precipitation_Fe3PO42', - reaction={'S_Fe2': -3*Fe_mw/Fe3PO42_mw/i_mass_Fe2, - 'S_IP': '?', - 'X_Fe3PO42': 1}, - ref_component='X_Fe3PO42', - conserved_for=('P',)) - precipitation.append(new_p) - - AlPO4_mw = cmps.X_AlPO4.chem_MW - Al_mw = cmps.S_Al.chem_MW - new_p = Process('precipitation_AlPO4', - reaction={'S_Al': -Al_mw/AlPO4_mw, - 'S_IP': '?', - 'X_AlPO4': 1}, - ref_component='X_AlPO4', - conserved_for=('P',)) - precipitation.append(new_p) - - self.extend(precipitation) - - gas_transfer = [] - for i in cls._biogas_IDs: - new_p = Process('%s_transfer' % i.lstrip('S_'), - reaction={i:-1}, - ref_component=i, - conserved_for=(), - parameters=()) - gas_transfer.append(new_p) - self.extend(gas_transfer) - self.compile(to_class=cls) - - stoichio_vals = (f_fa_li, f_bu_su, f_pro_su, f_ac_su, 1-f_bu_su-f_pro_su-f_ac_su, - f_va_aa, f_bu_aa, f_pro_aa, f_ac_aa, 1-f_va_aa-f_bu_aa-f_pro_aa-f_ac_aa, - f_ac_fa, 1-f_ac_fa, f_pro_va, f_ac_va, 1-f_pro_va-f_ac_va, - f_ac_bu, 1-f_ac_bu, f_ac_pro, 1-f_ac_pro, - Y_su, Y_aa, Y_fa, Y_c4, Y_pro, Y_ac, Y_h2, - f_ch_xb, f_pr_xb, f_li_xb, f_xI_xb, round(1.0-f_ch_xb-f_pr_xb-f_li_xb-f_xI_xb, 4), - f_va_pha, f_bu_pha, f_pro_pha, 1-f_va_pha-f_bu_pha-f_pro_pha, - 1-f_ac_pro, 1-f_ac_bu, 1-f_pro_va-f_ac_va, - Y_PO4, Y_hSRB, Y_aSRB, Y_pSRB, Y_c4SRB, - cmps.X_PP.i_K, cmps.X_PP.i_Mg, - cmps.S_S0.chem_MW, cmps.S_IS.chem_MW, - cmps.S_S0.i_mass, i_mass_IS, i_mass_Fe2) - - pH_limits = np.array([pH_limits_aa, pH_limits_ac, pH_limits_h2, - pH_limits_h2_SRB, pH_limits_ac_SRB, pH_limits_aa_SRB]).T - - ks = np.array((q_ch_hyd, q_pr_hyd, q_li_hyd, - k_su, k_aa, k_fa, k_c4, k_c4, k_pro, k_ac, k_h2, - b_su, b_aa, b_fa, b_c4, b_pro, b_ac, b_h2, # original ADM1 - q_pha, q_pha, q_pha, q_pha, b_pao, b_pp, b_pha, # P extension - k_hSRB, b_hSRB, k_aSRB, b_aSRB, k_pSRB, b_pSRB, k_c4SRB, k_c4SRB, b_c4SRB, # S extension - k_Fe3t2_h2, k_Fe3t2_h2, k_Fe3t2_is, k_Fe3t2_is, # Fe extension - q_aging_H, q_aging_L, q_Pcoprec, q_Pbinding, # HFO module - q_aging_H, q_aging_L, q_diss_H, q_diss_L)) - - Ks = np.array((K_su, K_aa, K_fa, K_c4, K_c4, K_pro, K_ac, K_h2, # original ADM1 - K_A, # P extension - K_hSRB, K_aSRB, K_pSRB, K_c4SRB, # S extension - K_Pbind, K_Pdiss)) # HFO module - K_so4 = np.array((K_so4_hSRB, K_so4_aSRB, K_so4_pSRB, K_so4_c4SRB)) - - KIs_h2 = np.array((KI_h2_fa, KI_h2_c4, KI_h2_c4, KI_h2_pro)) - KIs_h2s = np.array((KI_h2s_c4, KI_h2s_c4, KI_h2s_pro, KI_h2s_ac, KI_h2s_h2, - KI_h2s_hSRB, KI_h2s_aSRB, KI_h2s_pSRB, KI_h2s_c4SRB, KI_h2s_c4SRB)) - K_H_base = np.array(cls._K_H_base) - K_H_dH = np.array(cls._K_H_dH) - Ka_base = np.array([10**(-pKa) for pKa in pKa_base]) - Ka_dH = np.array(Ka_dH) - k_cryst = np.array(k_cryst) * 24 # converted to d^(-1) - n_cryst = np.array(n_cryst) - Ksp_base = np.array([10**(-pK) for pK in cls.pKsp_base]) - Ksp_dH = np.array(cls.Ksp_dH) - # root = TempState() - dct = self.__dict__ - dct.update(kwargs) - - dct['_parameters'] = dict(zip(cls._stoichio_params, stoichio_vals)) - self.set_rate_function(rhos_madm1) - self.rate_function._params = dict(zip(cls._kinetic_params, - [ks, Ks, K_PP, K_so4, - pH_limits, KS_IN*N_mw, KS_IP*P_mw, - KI_nh3, KIs_h2, KIs_h2s, - Ka_base, Ka_dH, K_H_base, K_H_dH, kLa, - k_cryst, n_cryst, Ksp_base, Ksp_dH, - cls.T_base, self._components, - # root, - ])) - return self \ No newline at end of file From a7abd5d6cc25e6b3cafb78f71d8cc34a4ad99f8e Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Wed, 23 Oct 2024 10:55:32 -0700 Subject: [PATCH 10/12] Revert "Update _adm1_p_extension.tsv" This reverts commit 8c2998eba231c4535e291f614d5ef0bae4e02fe9. --- qsdsan/data/process_data/_adm1_p_extension.tsv | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qsdsan/data/process_data/_adm1_p_extension.tsv b/qsdsan/data/process_data/_adm1_p_extension.tsv index 5ac6f4d8..5b91bb98 100644 --- a/qsdsan/data/process_data/_adm1_p_extension.tsv +++ b/qsdsan/data/process_data/_adm1_p_extension.tsv @@ -1,4 +1,4 @@ - S_su S_aa S_fa S_va S_bu S_pro S_ac S_h2 S_ch4 S_IC S_IN S_IP S_I X_ch X_pr X_li X_su X_aa X_fa X_c4 X_pro X_ac X_h2 X_I X_PHA X_PP X_PAO S_K S_Mg + S_su S_aa S_fa S_va S_bu S_pro S_ac S_h2 S_ch4 S_IC S_IN S_IP S_I X_ch X_pr X_li X_su X_aa X_fa X_c4 X_pro X_ac X_h2 X_I X_PHA X_PP X_PAO S_K S_Mg X_MeOH X_MeP hydrolysis_carbs 1 ? ? ? -1 hydrolysis_proteins 1 ? ? ? -1 hydrolysis_lipids 1-f_fa_li f_fa_li ? ? ? -1 From 7bf6ab28a2951688a643e41e1ab5b12f9412f288 Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Thu, 24 Oct 2024 16:28:05 -0700 Subject: [PATCH 11/12] added `PM2` doctest --- qsdsan/processes/_pm2.py | 35 ++++++++++++++++++++++++++++++++++- 1 file changed, 34 insertions(+), 1 deletion(-) diff --git a/qsdsan/processes/_pm2.py b/qsdsan/processes/_pm2.py index 24ebbafd..2d9f9709 100644 --- a/qsdsan/processes/_pm2.py +++ b/qsdsan/processes/_pm2.py @@ -610,6 +610,39 @@ class PM2(CompiledProcesses): 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]) + + >>> # Evaluate the rate of reaction at initial condition + >>> import numpy as np + >>> init_cond = { + ... 'X_CHL':2.81, + ... 'X_ALG':561.57, + ... 'X_CH':13.74, + ... 'X_LI':62.22, + ... 'S_CO2':30.0, + ... 'S_A':5.0, + ... 'S_F':5.0, + ... 'S_O2':20.36, + ... 'S_NH':25, + ... 'S_NO':9.30, + ... 'S_P':0.383, + ... 'X_N_ALG':3.62, + ... 'X_P_ALG':12.60, + ... } + >>> state_arr = np.append(cmps.kwarray(init_cond), [1000, 298, 112.6]) # flowrate, temperature, & irradiance + >>> pm2.rate_function(state_arr) + array([ 4.437, 142.045, 0.562, 0.562, 0.562, 2.48 , 304.319, + 193.505, 8.856, 24.737, 79.042, 2.093, 7.992, 67.419, + 186.095, 110.903, 5.076, 15.127, 32.669, 2.455, 9.375, + 45.795, 186.075, 110.903, 5.076, 15.126, 32.691, 2.456, + 9.378, 45.822]) + + >>> pm2.set_parameters(I_opt = 200) # Change optimal irradiance + >>> pm2.rate_function(state_arr) + array([ 4.437, 142.045, 0.562, 0.562, 0.562, 2.48 , 299.409, + 209.95 , 9.609, 34.175, 109.197, 2.093, 7.992, 67.419, + 171.898, 110.903, 5.076, 19.621, 42.373, 2.455, 9.375, + 45.795, 171.874, 110.903, 5.076, 19.618, 42.4 , 2.456, + 9.378, 45.822]) ''' _shared_params = ('Y_CH_PHO', 'Y_LI_PHO', 'Y_X_ALG_PHO', @@ -703,7 +736,7 @@ def set_parameters(self, **parameters): 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) + self.rate_function.set_params(**parameters) @property def Th_Q_N_min(self): From fee7d007318d99c8efa35c605293f417bb91e46b Mon Sep 17 00:00:00 2001 From: Joy Zhang Date: Tue, 29 Oct 2024 10:29:19 -0700 Subject: [PATCH 12/12] Update Systems.rst --- docs/source/Systems.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/Systems.rst b/docs/source/Systems.rst index 88a6c19e..d2de8daa 100644 --- a/docs/source/Systems.rst +++ b/docs/source/Systems.rst @@ -63,7 +63,7 @@ A variety of other sanitation and resource recovery systems have been developed #. Modular encapsulated two-stage anaerobic biological (METAB) system - * Manuscript: Zhang et al., Sustainable design of a modular anaerobic system for distributed energy recovery from industrial wastewaters, In Prep. + * Publication: `Zhang `_ et al., 2024. * `metab EXPOsan module `_ #. EcoRecover system: microalgae-based tertiary P recovery process