diff --git a/qsdsan/sanunits/_pfr.py b/qsdsan/sanunits/_pfr.py new file mode 100644 index 00000000..333dc601 --- /dev/null +++ b/qsdsan/sanunits/_pfr.py @@ -0,0 +1,80 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri May 24 08:18:43 2024 + +@author: joy_c +""" +import numpy as np +from numba import njit, cfunc + + +N = 6 +V = np.array([249, 2313, 2366, 2366, 4982, 4982]) +denom = np.diag(1/V) +f_in = np.array([[0, 0.8, 0.2, 0,0,0], + [1, 0, 0, 0, 0, 0]]) +internal_recycles = [(5, 2, 1.5e5)] # from, to, Q +Q_internal = np.zeros(N) +for i_from, i_to, qr in internal_recycles: + Q_internal[i_to: i_from] += qr + +DO = [0,0,0,0,2.,2.] +kLa = [0]*6 +DOsat = 8.0 + +#%% +from exposan import bsm1 +bsm1.load() +cmps = bsm1.PE.components +asm = bsm1.A1.suspended_growth_model +M_stoi = asm.stoichio_eval() +f_rho = asm.rate_function +DO_id = cmps.index('S_O') + +#%% +@njit +def dydt(t, QC_ins, QC, dQC_ins=None): + y = QC.reshape((N, len(QC)/N)) + Cs = y[:,:-1] + if any(DO): + for i in range(N): + if DO[i] > 0: Cs[i, DO_id] = DO[i] + elif any(kLa): + do = Cs[:, DO_id] + aer = kLa*(DOsat-do) + aer[aer < 0] = 0. + Qs = np.dot(QC_ins[:,-1], f_in).cumsum() + Q_internal + M_ins = f_in.T @ np.diag(QC_ins[:,-1]) @ QC_ins[:,:-1] # N * n_cmps + for i_from, i_to, qr in internal_recycles: + M_ins[i_to] += Cs[i_from] * qr + # rhos = np.apply_along_axis(f_rho, 1, Cs) # n_zone * n_process + rhos = np.vstack([f_rho(c) for c in Cs]) + rxn = rhos @ M_stoi + dy = np.zeros_like(y) + dy[:,:-1] = denom @ (M_ins - np.diag(Qs) @ Cs) + rxn + if any(DO): dy[:,DO_id] = 0. + elif any(kLa): dy[:,DO_id] += aer + return dy.flatten() + +#%% +import qsdsan.sanunits as su +from qsdsan import System +s = bsm1.sys.flowsheet.stream +u = bsm1.sys.flowsheet.unit +inf = s.wastewater +inf.unlink() +ras = s.RAS +ras.unlink() +s.treated.unlink() +AS = su.PFR('AS', ins=[inf, ras], outs=0-u.C1, V_tanks=[1000]*2+[1333]*3, + influent_fractions=[[1.0, 0,0,0,0], [1.0, 0,0,0,0]], + internal_recycles=[(4,0,55338)], kLa=[0,0,240,240,84], DO_ID='S_O', + suspended_growth_model=u.A1._model) +AS.set_init_conc( + S_I=30, S_S=5, X_I=1000, X_S=100, X_BH=500, X_BA=100, X_P=100, S_O=2, S_NO=20, + S_NH=2, S_ND=1, X_ND=1, S_ALK=84 + ) +sys = System('sys', path=(AS, u.C1)) + +#%% +sys.simulate(t_span=(0,50), method='BDF') diff --git a/qsdsan/sanunits/_suspended_growth_bioreactor.py b/qsdsan/sanunits/_suspended_growth_bioreactor.py index 955807b6..38d1aa18 100644 --- a/qsdsan/sanunits/_suspended_growth_bioreactor.py +++ b/qsdsan/sanunits/_suspended_growth_bioreactor.py @@ -19,11 +19,12 @@ 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', + 'PFR', ) def _add_aeration_to_growth_model(aer, model): @@ -732,16 +733,302 @@ def J_func(t, y): return J_func(*y) return (dC_dt, J_func) -# class PFR(SanUnit): +#%% +class PFR(SanUnit): + + _N_ins = 1 + _N_outs = 1 + _ins_size_is_fixed = False + _outs_size_is_fixed = True + + def __init__(self, ID='', ins=None, outs=(), thermo=None, init_with='WasteStream', + N_tanks_in_series=5, V_tanks=[1500, 1500, 3000, 3000, 3000], + influent_fractions=[[1.0, 0,0,0,0]], internal_recycles=[(4,0,35000),], + DO_setpoints=[], kLa=[0, 0, 120, 120, 60], DO_sat=8.0, + DO_ID='S_O2', suspended_growth_model=None, + isdynamic=True, exogenous_vars=(), **kwargs): + if exogenous_vars: + warn(f'currently exogenous dynamic variables are not supported in process simulation for {self.ID}') + SanUnit.__init__(self, ID, ins, outs, thermo, init_with, isdynamic=isdynamic, + exogenous_vars=exogenous_vars, **kwargs) + self.N_tanks_in_series = N_tanks_in_series + self.V_tanks = V_tanks + self.influent_fractions = influent_fractions + self.internal_recycles = internal_recycles + self.DO_setpoints = DO_setpoints + self.kLa = kLa + self.DO_sat = DO_sat + self.DO_ID = DO_ID + self.suspended_growth_model = suspended_growth_model + self._concs = None + self._Qs = self.V_tanks * 0 + + @property + def V_tanks(self): + '''[iterable[float]] Volumes of CSTRs in series [m3]''' + return self._Vs + @V_tanks.setter + def V_tanks(self, Vs): + if not iter(Vs): + raise TypeError(f'V_tanks must be an iterable, not {type(Vs).__name__}') + elif len(Vs) != self.N_tanks_in_series: + raise RuntimeError(f'cannot set the volumes of {self.N_tanks_in_series} tanks' + f'in series with {len(Vs)} value(s).') + else: self._Vs = np.asarray(Vs) + + @property + def influent_fractions(self): + '''[iterable[float]] Fractions of influents fed into different zones [unitless]''' + return self._f_ins + @influent_fractions.setter + def influent_fractions(self, fs): + fs = np.asarray(fs) + if len(fs.shape) == 1: + if len(self.ins) != 1: + raise RuntimeError(f'influent fractions must be a 2d array, with {len(self.ins)} row(s)' + f'and {self.N_tanks_in_series} columns') + if fs.shape[0] != self.N_tanks_in_series: + raise RuntimeError('cannot set the fractions of influent fed into ' + f'{self.N_tanks_in_series} tanks in series with {len(fs)} value(s).') + fs = fs.reshape(1, len(fs)) + elif len(fs.shape) > 2: + raise RuntimeError(f'influent fractions must be a 2d array, with {len(self.ins)} row(s)' + f'and {self.N_tanks_in_series} columns') + if (fs < 0).any(): + raise ValueError('influent fractions must not have negative value(s)') + for row in fs: + rowsum = row.sum() + if rowsum != 1: row /= rowsum + self._f_ins = fs + + @property + def internal_recycles(self): + '''[list[3-tuple[int, int, float]]] A list of three-element tuples (i, j, Q) indicating internal recycling + streams from zone i to zone j with a flowrate of Q [m3/d], respectively. + Indices i,j start from 0 not 1.''' + return self._rcy + @internal_recycles.setter + def internal_recycles(self, rcy): + isa = isinstance + if isa(rcy, tuple): + if len(rcy) != 3: + raise TypeError('internal recycles must be indicated by a list of 3-tuples') + rcy = [rcy] + elif isa(rcy, list): + for row in rcy: + if len(row) != 3: + raise TypeError('internal recycles must be indicated by a list of 3-tuples') + else: + raise TypeError('internal recycles must be indicated by a list of 3-tuples') + _rcy = [] + for row in rcy: + i, j, q = row + _rcy.append((int(i), int(j), q)) + self._rcy = _rcy + + @property + def DO_setpoints(self): + '''[iterable[float]] Dissolve oxygen setpoints of CSTRs in series [mg-O2/L]. + 0 is treated as no active aeration. DO setpoints take priority over kLa values.''' + return self._DOs + @DO_setpoints.setter + def DO_setpoints(self, DOs): + if not iter(DOs): + raise TypeError(f'DO setpoints must be an iterable, not {type(DOs).__name__}') + elif 0 < len(DOs) != self.N_tanks_in_series: + raise RuntimeError(f'cannot set the DO setpoints of {self.N_tanks_in_series} tanks' + f'in series with {len(DOs)} value(s).') + else: self._DOs = np.asarray(DOs) + + @property + def kLa(self): + '''[iterable[float]] Aeration kLa values of CSTRs in series [d^(-1)]. If DO + setpoints are specified, kLa values would be ignored in process simulation.''' + return self._kLa + @kLa.setter + def kLa(self, ks): + if not iter(ks): + raise TypeError(f'V_tanks must be an iterable, not {type(ks).__name__}') + elif len(ks) != self.N_tanks_in_series: + raise RuntimeError(f'cannot set kLa of {self.N_tanks_in_series} tanks' + f'in series with {len(ks)} value(s).') + else: self._kLa = np.asarray(ks) + + @property + def suspended_growth_model(self): + '''[:class:`CompiledProcesses` or NoneType] Suspended growth model.''' + return self._model + + @suspended_growth_model.setter + def suspended_growth_model(self, model): + if isinstance(model, CompiledProcesses) or model is None: self._model = model + else: raise TypeError(f'suspended_growth_model must be one of the following ' + f'types: CompiledProesses, NoneType. Not {type(model)}') + + @property + def DO_ID(self): + '''[str] The `Component` ID for dissolved oxygen used in the suspended growth model and the aeration model.''' + return self._DO_ID + + @DO_ID.setter + def DO_ID(self, doid): + if doid not in self.components.IDs: + raise ValueError(f'DO_ID must be in the set of `CompiledComponents` used to set thermo, ' + f'i.e., one of {self.components.IDs}.') + self._DO_ID = doid + + + def _run(self): + out, = self.outs + out.mix_from(self.ins) + + @property + def state(self): + '''The state of the PFR, including component concentrations [mg/L] and flow rate [m^3/d] for each zone.''' + if self._state is None: return None + else: + N = self.N_tanks_in_series + y = self._state.copy() + y = y.reshape((N, int(len(y)/N))) + y = pd.DataFrame(y, columns=self._state_header) + return y + + def set_init_conc(self, concentrations=None, i_zone=None, **kwargs): + '''set the initial concentrations [mg/L] of specific zones.''' + isa = isinstance + cmps = self.components + N = self.N_tanks_in_series + if self._concs is None: self._concs = np.zeros((N, len(cmps))) + if concentrations is None: + concs = cmps.kwarray(kwargs) + if i_zone is None: + self._concs[:] = concs + else: + self._concs[i_zone] = concs + elif isa(concentrations, pd.DataFrame): + concentrations.index = range(N) + dct = concentrations.to_dict('index') + for i, concs in dct: + self._concs[i] = cmps.kwarray(concs) + elif isa(concentrations, np.ndarray): + if concentrations.shape != self._concs.shape: + raise RuntimeError(f'cannot set the concentrations of {len(cmps)} ' + f'components across {N} with a {concentrations.shape} array') + self._concs = concentrations + else: + raise TypeError('specify initial concentrations with pandas.DataFrame, numpy.ndarray' + 'or " = value" kwargs') + + def _init_state(self): + self._state_header = ny = list(self.components.IDs) + ['Q'] + self._Qs_idx = list(len(ny) * np.arange(1, 1+self.N_tanks_in_series) - 1) + out, = self.outs + y = np.zeros((self.N_tanks_in_series, len(ny))) + if self._concs is not None: + y[:,:-1] = self._concs + else: + y[:,:-1] = out.conc + y[:,-1] = out.F_vol*24 + self._state = y.flatten() + self._dstate = self._state * 0. + + def _update_state(self): + out, = self.outs + n_cmp = len(self.components) + self._state[self._Qs_idx] = self._Qs + out.state[:-1] = self._state[-(n_cmp+1):-1] + out.state[-1] = sum(ws.state[-1] for ws in self.ins) + + def _update_dstate(self): + out, = self.outs + n_cmp = len(self.components) + out.dstate[:-1] = self._dstate[-(n_cmp+1):-1] + out.dstate[-1] = sum(ws.dstate[-1] for ws in self.ins) -# _N_ins = 1 -# _N_outs = 2 + @property + def ODE(self): + if self._ODE is None: + self._compile_ODE() + return self._ODE -# def __init__(self, ID='', ins=None, outs=(), **kwargs): -# SanUnit.__init__(self, ID, ins, outs) + def _compile_ODE(self): + _Qs = self._Qs + _dstate = self._dstate + _update_dstate = self._update_dstate + N = self.N_tanks_in_series + _1_ov_V = np.diag(1/self.V_tanks) + f_in = self.influent_fractions + DO = self.DO_setpoints + kLa = self.kLa + rcy = self.internal_recycles + DO_idx = self.components.index(self.DO_ID) + DOsat = self.DO_sat + Q_internal = np.zeros(N) + for i_from, i_to, qr in rcy: + Q_internal[i_to: i_from] += qr -# def _run(self, steady_state=True): -# pass + if self._model is None: + warn(f'{self.ID} was initialized without a suspended growth model, ' + f'and thus run as a non-reactive unit') + Rs = lambda Cs: 0. + else: + f_rho = self._model.rate_function + M_stoi = self._model.stoichio_eval() + # @njit + def Rs(Cs): + rhos = np.vstack([f_rho(c) for c in Cs]) + rxn = rhos @ M_stoi + return rxn + + if any(DO): + aerated_zones = (DO > 0) + aerated_DO = DO[aerated_zones] + # @njit + def dy_dt(t, QC_ins, QC, dQC_ins): + y = QC.reshape((N, int(len(QC)/N))) + Cs = y[:,:-1] + Cs[aerated_zones, DO_idx] = aerated_DO + _Qs[:] = Qs = np.dot(QC_ins[:,-1], f_in).cumsum() + Q_internal + M_ins = f_in.T @ np.diag(QC_ins[:,-1]) @ QC_ins[:,:-1] # N * n_cmps + M_outs = np.diag(Qs) @ Cs + M_ins[1:] += M_outs[:-1] + for i_from, i_to, qr in rcy: + M_rcy = Cs[i_from] * qr + M_ins[i_to] += M_rcy + if i_from + 1 < N: M_ins[i_from+1] -= M_rcy + rxn = Rs(Cs) + dy = np.zeros_like(y) + dy[:,:-1] = _1_ov_V @ (M_ins - M_outs) + rxn + dy[aerated_zones, DO_idx] = 0. + _dstate[:] = dy.flatten() + _update_dstate() + + else: + if not any(kLa): kLa = np.zeros(N) + # @njit + def dy_dt(t, QC_ins, QC, dQC_ins): + # breakpoint() + y = QC.reshape((N, int(len(QC)/N))) + Cs = y[:,:-1] + do = Cs[:, DO_idx] + aer = kLa*(DOsat-do) + # aer[aer < 0] = 0. + _Qs[:] = Qs = np.dot(QC_ins[:,-1], f_in).cumsum() + Q_internal + M_ins = f_in.T @ np.diag(QC_ins[:,-1]) @ QC_ins[:,:-1] # N * n_cmps + M_outs = np.diag(Qs) @ Cs + M_ins[1:] += M_outs[:-1] + for i_from, i_to, qr in rcy: + M_rcy = Cs[i_from] * qr + M_ins[i_to] += M_rcy + if i_from + 1 < N: M_ins[i_from+1] -= M_rcy + rxn = Rs(Cs) + dy = np.zeros_like(y) + dy[:,:-1] = _1_ov_V @ (M_ins - M_outs) + rxn + dy[:,DO_idx] += aer + _dstate[:] = dy.flatten() + _update_dstate() + + self._ODE = dy_dt -# def _design(self): -# pass \ No newline at end of file + def _design(self): + pass \ No newline at end of file