Skip to content

Commit

Permalink
Merge branch 'metro-bsm2' into metro
Browse files Browse the repository at this point in the history
  • Loading branch information
@RaiSaumitra
RaiSaumitra committed May 10, 2024
2 parents d5b26ba + 61cd5ec commit e05a275
Show file tree
Hide file tree
Showing 10 changed files with 404 additions and 286 deletions.
19 changes: 14 additions & 5 deletions qsdsan/_process.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1192,8 +1192,12 @@ def stoichiometry(self):
if isa(stoichio, list) and len(v_params) > 0:
stoichio_vals = []
for row in stoichio:
stoichio_vals.append([v.subs(v_params) if not isa(v, (float, int)) else v for v in row])
try: stoichio_vals = np.asarray(stoichio_vals, dtype=float)
# stoichio_vals.append([v.subs(v_params) if not isa(v, (float, int)) else v for v in row])
stoichio_vals.append([v.evalf(subs=v_params) if not isa(v, (float, int)) else v for v in row])
try:
stoichio_vals = np.asarray(stoichio_vals, dtype=float)
#!!! round-off error
# stoichio_vals[np.abs(stoichio_vals) < 2.22044604925e-16] = 0.0
except TypeError: pass
return pd.DataFrame(stoichio_vals, index=self.IDs, columns=self._components.IDs)
else: return pd.DataFrame(stoichio, index=self.IDs, columns=self._components.IDs)
Expand All @@ -1202,11 +1206,16 @@ def _lambdify_stoichio(self):
dct = self._dyn_params
dct_vals = self._parameters
if dct:
sbs = [i.symbol for i in dct.values()]
lamb = lambdify(sbs, self._stoichiometry, 'numpy')
static_params = {k:v for k,v in dct_vals.items() if k not in dct}
stoichio = []
isa = isinstance
for row in self._stoichiometry:
stoichio.append([v.evalf(subs=static_params) if not isa(v, (float, int)) else v for v in row])
sbs = [symbols(i) for i in dct.keys()]
lamb = lambdify(sbs, stoichio, 'numpy')
arr = np.empty((self.size, len(self._components)))
def f():
v = [v for k,v in dct_vals.items() if k in dct.keys()]
v = [dct_vals[k] for k in dct.keys()]
arr[:,:] = lamb(*v)
return arr
self.__dict__['_stoichio_lambdified'] = f
Expand Down
3 changes: 2 additions & 1 deletion qsdsan/_sanunit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -269,7 +269,8 @@ def _convert_stream(self, strm_inputs, streams, init_with, ins_or_outs):
elif v == 'ss':
converted.append(SanStream.from_stream(stream=s))
else:
converted.append(WasteStream.from_stream(stream=s))
if isa(s, WasteStream): converted.append(s)
else: converted.append(WasteStream.from_stream(stream=s))

diff = len(converted) + len(missing) - len(streams)
if diff != 0:
Expand Down
19 changes: 12 additions & 7 deletions qsdsan/_waste_stream.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -452,15 +452,16 @@ def _wastestream_info(self, details=True, concentrations=None, N=15):
_ws_info += int(bool(self.pH))*f' pH : {self.pH:.1f}\n'
_ws_info += int(bool(self.SAlk))*f' Alkalinity : {self.SAlk:.1f} mg/L\n'
if details:
_ws_info += int(bool(self.COD)) *f' COD : {self.COD:.1f} mg/L\n'
_ws_info += int(bool(self.BOD)) *f' BOD : {self.BOD:.1f} mg/L\n'
_ws_info += int(bool(self.TC)) *f' TC : {self.TC:.1f} mg/L\n'
_ws_info += int(bool(self.TOC)) *f' TOC : {self.TOC:.1f} mg/L\n'
_ws_info += int(bool(self.TN)) *f' TN : {self.TN:.1f} mg/L\n'
_ws_info += int(bool(self.COD)) *f' COD : {self.COD:.1f} mg/L\n'
_ws_info += int(bool(self.BOD)) *f' BOD : {self.BOD:.1f} mg/L\n'
_ws_info += int(bool(self.TC)) *f' TC : {self.TC:.1f} mg/L\n'
_ws_info += int(bool(self.TOC)) *f' TOC : {self.TOC:.1f} mg/L\n'
_ws_info += int(bool(self.TN)) *f' TN : {self.TN:.1f} mg/L\n'
# `TKN` not included as the users need to define that to include in TKN calculation
# _ws_info += int(bool(self.TKN)) *f' TKN : {self.TKN:.1f} mg/L\n'
_ws_info += int(bool(self.TP)) *f' TP : {self.TP:.1f} mg/L\n'
_ws_info += int(bool(self.TK)) *f' TK : {self.TK:.1f} mg/L\n'
_ws_info += int(bool(self.TP)) *f' TP : {self.TP:.1f} mg/L\n'
_ws_info += int(bool(self.TK)) *f' TK : {self.TK:.1f} mg/L\n'
_ws_info += int(bool(self.get_TSS())) *f' TSS : {self.get_TSS():.1f} mg/L\n'
# _ws_info += int(bool(self.charge))*f' charge : {self.charge:.1f} mmol/L\n'
else:
_ws_info += ' ...\n'
Expand Down Expand Up @@ -729,6 +730,10 @@ def _liq_sol_properties(self, prop, value):
def pH(self):
'''[float] pH, unitless.'''
return self._liq_sol_properties('pH', 7.)
@pH.setter
def pH(self, ph):
if self.phase != 'g':
self._pH = ph

@property
def SAlk(self):
Expand Down
133 changes: 88 additions & 45 deletions qsdsan/processes/_adm1.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,9 @@
'non_compet_inhibit', 'substr_inhibit',
'T_correction_factor',
'pH_inhibit', 'Hill_inhibit',
'rhos_adm1')
'rhos_adm1',
'solve_pH',
'dydt_Sh2_AD', 'grad_dydt_Sh2_AD')

_path = ospath.join(data_path, 'process_data/_adm1.tsv')
_load_components = settings.get_default_chemicals
Expand All @@ -36,8 +38,10 @@
# ADM1-specific components
# =============================================================================

C_mw = get_mw({'C':1})
N_mw = get_mw({'N':1})
# C_mw = get_mw({'C':1})
# N_mw = get_mw({'N':1})
C_mw = 12
N_mw = 14

def create_adm1_cmps(set_thermo=True):
cmps_all = Components.load_default()
Expand Down Expand Up @@ -179,25 +183,25 @@ def create_adm1_cmps(set_thermo=True):
def non_compet_inhibit(Si, Ki):
return Ki/(Ki+Si)

def grad_non_compet_inhibit(Si, Ki):
return -Ki/(Ki+Si)**2

def substr_inhibit(Si, Ki):
return Si/(Ki+Si)

def grad_substr_inhibit(Si, Ki):
return Ki/(Ki+Si)**2

def mass2mol_conversion(cmps):
'''conversion factor from kg[measured_as]/m3 to mol[component]/L'''
return cmps.i_mass / cmps.chem_MW

# def T_correction_factor(T1, T2, theta):
# return np.exp(theta * (T2-T1))

def T_correction_factor(T1, T2, delta_H):
"""compute temperature correction factor for equilibrium constants based on
the Van't Holf equation."""
if T1 == T2: return 1
return np.exp(delta_H/(R*100) * (1/T1 - 1/T2)) # R converted to SI

# def calc_Kas(pKas, T_base, T_op, theta):
# pKas = np.asarray(pKas)
# return 10**(-pKas) * T_correction_factor(T_base, T_op, theta)

def acid_base_rxn(h_ion, weak_acids_tot, Kas):
# h, nh4, hco3, ac, pr, bu, va = mols
Expand All @@ -218,8 +222,6 @@ def fprime_abr(h_ion, weak_acids_tot, Kas):

def pH_inhibit(pH, ul, ll, lower_only=True):
if lower_only:
# if pH >= ul: return 1
# else: return exp(-3 * ((pH-ul)/(ul-ll))**2)
low_by = np.minimum(pH-ul, 0)
return np.exp(-3 * (low_by/(ul-ll))**2)
else:
Expand All @@ -233,11 +235,22 @@ def Hill_inhibit(H_ion, ul, ll):
rhos = np.zeros(22) # 22 kinetic processes
Cs = np.empty(19)

def rhos_adm1(state_arr, params):
def solve_pH(state_arr, Ka, unit_conversion):
cmps_in_M = state_arr[:27] * unit_conversion
weak_acids = cmps_in_M[[24, 25, 10, 9, 6, 5, 4, 3]]
h = brenth(acid_base_rxn, 1e-14, 1.0,
args=(weak_acids, Ka),
xtol=1e-12, maxiter=100)
nh3 = Ka[1] * weak_acids[2] / (Ka[1] + h)
co2 = weak_acids[3] - Ka[2] * weak_acids[3] / (Ka[2] + h)
return h, nh3, co2

rhos_adm1 = lambda state_arr, params: _rhos_adm1(state_arr, params, h=None)

def _rhos_adm1(state_arr, params, h=None):
ks = params['rate_constants']
Ks = params['half_sat_coeffs']
cmps = params['components']
# n = len(cmps)
pH_ULs = params['pH_ULs']
pH_LLs = params['pH_LLs']
KS_IN = params['KS_IN']
Expand All @@ -250,26 +263,20 @@ def rhos_adm1(state_arr, params):
kLa = params['kLa']
T_base = params['T_base']
root = params['root']
if 'unit_conv' in params:
unit_conversion = params['unit_conv']
else:
unit_conversion = params['unit_conv'] = mass2mol_conversion(cmps)

# Cs_ids = cmps.indices(['X_c', 'X_ch', 'X_pr', 'X_li', 'X_su', 'X_aa',
# 'X_fa', 'X_c4', 'X_c4', 'X_pro', 'X_ac', 'X_h2',
# 'X_su', 'X_aa', 'X_fa', 'X_c4', 'X_pro', 'X_ac', 'X_h2'])
# Cs = state_arr[Cs_ids]
Cs[:8] = state_arr[12:20]
Cs[8:12] = state_arr[19:23]
Cs[12:] = state_arr[16:23]
# substrates_ids = cmps.indices(['S_su', 'S_aa', 'S_fa', 'S_va',
# 'S_bu', 'S_pro', 'S_ac', 'S_h2'])
# substrates = state_arr[substrates_ids]

substrates = state_arr[:8]
# S_va, S_bu, S_h2, S_IN = state_arr[cmps.indices(['S_va', 'S_bu', 'S_h2', 'S_IN'])]
# S_va, S_bu, S_h2, S_ch4, S_IC, S_IN = state_arr[[3,4,7,8,9,10]]
S_va, S_bu, S_h2, S_IN = state_arr[[3,4,7,10]]
unit_conversion = mass2mol_conversion(cmps)
cmps_in_M = state_arr[:27] * unit_conversion
weak_acids = cmps_in_M[[24, 25, 10, 9, 6, 5, 4, 3]]

T_op = state_arr[-1]
# Ka, KH = T_corrected_params(T_op, params)
if T_op == T_base:
Ka = Kab
KH = KHb / unit_conversion[7:10]
Expand All @@ -286,22 +293,19 @@ def rhos_adm1(state_arr, params):

biogas_S = state_arr[7:10].copy()
biogas_p = R * T_op * state_arr[27:30]
# Kas = Kab * T_correction_factor(T_base, T_op, Ka_dH)
# KH = KHb * T_correction_factor(T_base, T_op, KH_dH) / unit_conversion[7:10]

rhos[:-3] = ks * Cs
Monod = substr_inhibit(substrates, Ks)
rhos[4:12] *= Monod
if S_va > 0: rhos[7] *= 1/(1+S_bu/S_va)
if S_bu > 0: rhos[8] *= 1/(1+S_va/S_bu)

h = brenth(acid_base_rxn, 1e-14, 1.0,
args=(weak_acids, Ka),
xtol=1e-12, maxiter=100)
# h = 10**(-7.46)

nh3 = Ka[1] * weak_acids[2] / (Ka[1] + h)
co2 = weak_acids[3] - Ka[2] * weak_acids[3] / (Ka[2] + h)
if h is None:
h, nh3, co2 = solve_pH(state_arr, Ka, unit_conversion)
else:
nh3 = Ka[1] * S_IN * unit_conversion[10] / (Ka[1] + h)
S_IC = state_arr[9] * unit_conversion[9]
co2 = S_IC - Ka[2] * S_IC / (Ka[2] + h)
biogas_S[-1] = co2 / unit_conversion[9]

Iph = Hill_inhibit(h, pH_ULs, pH_LLs)
Expand All @@ -310,8 +314,6 @@ def rhos_adm1(state_arr, params):
Inh3 = non_compet_inhibit(nh3, KI_nh3)
rhos[4:12] *= Iph * Iin
rhos[6:10] *= Ih2
# rhos[4:12] *= Hill_inhibit(h, pH_ULs, pH_LLs) * substr_inhibit(S_IN, KS_IN)
# rhos[6:10] *= non_compet_inhibit(S_h2, KIs_h2)
rhos[10] *= Inh3
root.data = {
'pH':-np.log10(h),
Expand All @@ -323,9 +325,47 @@ def rhos_adm1(state_arr, params):
'rhos':rhos[4:12].copy()
}
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):
state_arr[7] = S_h2
Q = state_arr[30]
rxn = _rhos_adm1(state_arr, params, h=h)
stoichio = f_stoichio(state_arr) # should return the stoichiometric coefficients of S_h2 for all processes
return Q/V_liq*(S_h2_in - S_h2) + np.dot(rxn, stoichio)

grad_rhos = np.zeros(5)
X_bio = np.zeros(5)
def grad_dydt_Sh2_AD(S_h2, state_arr, h, params, f_stoichio, V_liq, S_h2_in):
state_arr[7] = S_h2
ks = params['rate_constants'][[6,7,8,9,11]]
Ks = params['half_sat_coeffs'][2:6]
K_h2 = params['half_sat_coeffs'][7]
pH_ULs = params['pH_ULs']
pH_LLs = params['pH_LLs']
KS_IN = params['KS_IN']
KIs_h2 = params['KIs_h2']
kLa = params['kLa']

X_bio[:] = state_arr[[18,19,19,20,22]]
substrates = state_arr[2:6]
S_va, S_bu, S_IN = state_arr[[3,4,10]]
Iph = Hill_inhibit(h, pH_ULs, pH_LLs)[[2,3,4,5,7]]
Iin = substr_inhibit(S_IN, KS_IN)
grad_Ih2 = grad_non_compet_inhibit(S_h2, KIs_h2)

grad_rhos[:] = ks * X_bio * Iph * Iin
grad_rhos[:-1] *= substr_inhibit(substrates, Ks) * grad_Ih2
if S_va > 0: rhos[1] *= 1/(1+S_bu/S_va)
if S_bu > 0: rhos[2] *= 1/(1+S_va/S_bu)

grad_rhos[-1] *= grad_substr_inhibit(S_h2, K_h2)
stoichio = f_stoichio(state_arr)

Q = state_arr[30]
return -Q/V_liq + np.dot(grad_rhos, stoichio[[6,7,8,9,11]]) + kLa*stoichio[-3]


#%%
# =============================================================================
# ADM1 class
Expand Down Expand Up @@ -569,9 +609,12 @@ def __new__(cls, components=None, path=None, N_xc=2.686e-3, N_I=4.286e-3, N_aa=7
**kwargs):

cmps = _load_components(components)
cmps.X_c.i_N = N_xc * N_mw
cmps.X_I.i_N = cmps.S_I.i_N = N_I * N_mw
cmps.S_aa.i_N = cmps.X_pr.i_N = N_aa * N_mw
cmps.X_c.i_N = round(N_xc * N_mw, 4)
cmps.X_I.i_N = cmps.S_I.i_N = round(N_I * N_mw, 4)
cmps.S_aa.i_N = cmps.X_pr.i_N = round(N_aa * N_mw, 4)
# cmps.X_c.i_N = N_xc * N_mw
# cmps.X_I.i_N = cmps.S_I.i_N = N_I * N_mw
# cmps.S_aa.i_N = cmps.X_pr.i_N = N_aa * N_mw

if not path: path = _path
self = Processes.load_from_file(path,
Expand All @@ -591,11 +634,11 @@ def __new__(cls, components=None, path=None, N_xc=2.686e-3, N_I=4.286e-3, N_aa=7
self.extend(gas_transfer)
self.compile(to_class=cls)

stoichio_vals = (f_ch_xc, f_pr_xc, f_li_xc, f_xI_xc, 1-f_ch_xc-f_pr_xc-f_li_xc-f_xI_xc,
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,
stoichio_vals = (f_ch_xc, f_pr_xc, f_li_xc, f_xI_xc, 1.0-f_ch_xc-f_pr_xc-f_li_xc-f_xI_xc,
f_fa_li, f_bu_su, f_pro_su, f_ac_su, 1.0-f_bu_su-f_pro_su-f_ac_su,
f_va_aa, f_bu_aa, f_pro_aa, f_ac_aa, 1.0-f_va_aa-f_bu_aa-f_pro_aa-f_ac_aa,
f_ac_fa, 1.0-f_ac_fa, f_pro_va, f_ac_va, 1.0-f_pro_va-f_ac_va,
f_ac_bu, 1.0-f_ac_bu, f_ac_pro, 1.0-f_ac_pro,
Y_su, Y_aa, Y_fa, Y_c4, Y_pro, Y_ac, Y_h2)
pH_LLs = np.array([pH_limits_aa[0]]*6 + [pH_limits_ac[0], pH_limits_h2[0]])
pH_ULs = np.array([pH_limits_aa[1]]*6 + [pH_limits_ac[1], pH_limits_h2[1]])
Expand Down
11 changes: 10 additions & 1 deletion qsdsan/sanunits/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,16 @@
Please refer to https://github.com/QSD-Group/QSDsan/blob/main/LICENSE.txt
for license details.
'''

# %%
from numba import njit
@njit(cache=True)
def dydt_cstr(QC_ins, QC, V, _dstate):
Q_ins = QC_ins[:, -1]
C_ins = QC_ins[:, :-1]
_dstate[-1] = 0
_dstate[:-1] = (Q_ins @ C_ins - sum(Q_ins)*QC[:-1])/V

#%%
# **NOTE** PLEASE ORDER THE MODULES ALPHABETICALLY #

# Units that do not rely on other units
Expand Down
Loading

0 comments on commit e05a275

Please sign in to comment.