Changes to allow customized properties and impact item creation

qsdsan/_component.py

@@ -233,6 +233,10 @@ def __new__(cls, ID, search_ID=None, formula=None, phase=None, measured_as=None,
                 self._chem_MW = molecular_weight(self.atoms)
         if phase: lock_phase(self, phase)
 
+        self._particle_size = particle_size
+        self._degradability = degradability
+        self._organic = organic
+        self.description = description
         self._measured_as = measured_as
         self.i_mass = i_mass
         self.i_C = i_C
@@ -245,10 +249,7 @@ def __new__(cls, ID, search_ID=None, formula=None, phase=None, measured_as=None,
         self.f_BOD5_COD = f_BOD5_COD
         self.f_uBOD_COD = f_uBOD_COD
         self.f_Vmass_Totmass = f_Vmass_Totmass
-        self._particle_size = particle_size
-        self._degradability = degradability
-        self._organic = organic
-        self.description = description
+
         if not self.MW and not self.formula: self.MW = 1.
         self.i_COD = i_COD
         self.i_NOD = i_NOD
@@ -340,7 +341,8 @@ def i_mass(self, i):
                     chem_charge = charge_from_formula(self.formula)
                     Cr2O7 = - cod_test_stoichiometry(self.atoms, chem_charge)['Cr2O7-2']
                     cod = Cr2O7 * 1.5 * molecular_weight({'O':2})
-                    i = self.chem_MW/cod
+                    try: i = self.chem_MW/cod
+                    except: breakpoint()
                 elif self.measured_as:
                     raise AttributeError(f'Must specify i_mass for component {self.ID} '
                                          f'measured as {self.measured_as}.')
@@ -532,7 +534,7 @@ def i_COD(self):
     def i_COD(self, i):
         if i is not None: self._i_COD = check_return_property('i_COD', i)
         else:
-            if self.organic or self.formula in ('H2', 'O2', 'N2', 'NO2-', 'NO3-'):
+            if self.organic or self.formula in ('H2', 'O2', 'N2', 'NO2-', 'NO3-', 'H2S', 'S'):
                 if self.measured_as == 'COD': self._i_COD = 1.
                 elif not self.atoms:
                     raise AttributeError(f"Must specify `i_COD` for organic component {self.ID}, "
@@ -792,6 +794,10 @@ def from_chemical(cls, ID, chemical=None, formula=None, phase=None, measured_as=
                            new.Tm, new.Tb, new.eos, new.phase_ref, new.S0)
         TDependentProperty.RAISE_PROPERTY_CALCULATION_ERROR = True
 
+        new.description = description
+        new._particle_size = particle_size
+        new._degradability = degradability
+        new._organic = organic
         new._measured_as = measured_as
         new.i_mass = i_mass
         new.i_C = i_C
@@ -804,10 +810,7 @@ def from_chemical(cls, ID, chemical=None, formula=None, phase=None, measured_as=
         new.f_BOD5_COD = f_BOD5_COD
         new.f_uBOD_COD = f_uBOD_COD
         new.f_Vmass_Totmass = f_Vmass_Totmass
-        new.description = description
-        new._particle_size = particle_size
-        new._degradability = degradability
-        new._organic = organic
+
         new.i_COD = i_COD
         new.i_NOD = i_NOD
         return new

qsdsan/_impact_item.py

@@ -383,8 +383,12 @@ def get_item(cls, ID):
     @classmethod
     def _load_from_df(cls, name, df):
         if name.lower() == 'info':
+            if 'kind' not in df.columns: df['kind'] = 'ImpactItem'
             for num in df.index:
-                new = cls.__new__(cls)
+                kind = df.iloc[num].kind.lower()
+                if kind == 'streamimpactitem':
+                    new = StreamImpactItem.__new__(StreamImpactItem)
+                else: new = cls.__new__(cls)
                 new.__init__(ID=df.iloc[num].ID,
                              functional_unit=df.iloc[num].functional_unit)
         else:
@@ -409,8 +413,9 @@ def load_from_file(cls, path_or_dict, index_col=None):
 
         This Excel should have multiple sheets:
 
-            - The "info" sheet should have two columns: "ID" (e.g., Cement) \
-            and "functional_unit" (e.g., kg) of different impact items.
+            - The "info" sheet should have three columns: "ID" (e.g., Cement) \
+             "functional_unit" (e.g., kg), and "kind" ("ImpactItem" or "StreamImpactItem") 
+             of different impact items.
 
             - The remaining sheets should contain characterization factors of \
             impact indicators.

qsdsan/_process.py

@@ -706,8 +706,9 @@ def _normalize_stoichiometry(self, new_ref):
             self._stoichiometry = [v/factor for v in stoich]
 
     def _normalize_rate_eq(self, new_ref):
-        factor = abs(self._stoichiometry[self._components.index(str(new_ref))])
-        self._rate_equation *= factor
+        if self._rate_equation:
+            factor = abs(self._stoichiometry[self._components.index(str(new_ref))])
+            self._rate_equation *= factor
 
     def show(self):
         info = f"Process: {self.ID}"

qsdsan/_waste_stream.py

@@ -40,12 +40,13 @@
 _defined_composite_vars = ('COD', 'BOD5', 'BOD', 'uBOD', 'NOD', 'ThOD', 'cnBOD',
                            'C', 'N', 'P', 'K', 'Mg', 'Ca', 'solids', 'charge')
 
-_common_composite_vars = ('_COD', '_BOD', '_uBOD', '_TC', '_TOC', '_TN',
+_common_composite_vars = ('_COD', '_BOD', '_uBOD', '_ThOD', '_cnBOD',
+                          '_TC', '_TOC', '_TN',
                           '_TKN', '_TP', '_TK', '_TMg', '_TCa',
-                          '_dry_mass', '_charge', '_ThOD', '_cnBOD')
+                          '_dry_mass', '_charge',)
 
 _ws_specific_slots = (*_common_composite_vars,
-                      '_pH', '_SAlk', '_ratios',
+                      '_pH', '_SAlk', '_ratios', 'additional_properties',
                       # '_stream_impact_item', (pls keep this here, might be useful in debugging)
                       '_state', '_dstate', '_scope')
 
@@ -249,7 +250,9 @@ class WasteStream(SanStream):
         TO BE IMPLEMENTED
     stream_impact_item : :class:`StreamImpactItem`
         The :class:`StreamImpactItem` this stream is linked to.
-    component_flows : kwargs
+    additional_properties : dict
+        Additional properties (e.g., turbidity).
+    component_flows : dict
         Component flow data.
 
 
@@ -271,20 +274,28 @@ def __init__(self, ID='', flow=(), phase='l', T=298.15, P=101325.,
                  pH=7., SAlk=2.5, COD=None, BOD=None, uBOD=None,
                  ThOD=None, cnBOD=None,
                  TC=None, TOC=None, TN=None, TKN=None, TP=None, TK=None,
-                 TMg=None, TCa=None, dry_mass=None, charge=None, ratios=None,
-                 stream_impact_item=None, **component_flows):
+                 TMg=None, TCa=None,
+                 dry_mass=None, charge=None, ratios=None,
+                 stream_impact_item=None, additional_properties={},
+                 **component_flows):
         SanStream.__init__(self=self, ID=ID, flow=flow, phase=phase, T=T, P=P,
                            units=units, price=price, thermo=thermo,
                            stream_impact_item=stream_impact_item, **component_flows)
 
-        self._init_ws(pH, SAlk, COD, BOD, uBOD, TC, TOC, TN, TKN,
-                      TP, TK, TMg, TCa, ThOD, cnBOD, dry_mass, charge, ratios)
+        self._init_ws(pH=pH, SAlk=SAlk, COD=COD, BOD=BOD, uBOD=uBOD,
+                      ThOD=ThOD, cnBOD=cnBOD, 
+                      TC=TC, TOC=TOC, TN=TN, TKN=TKN,
+                      TP=TP, TK=TK, TMg=TMg, TCa=TCa,
+                      dry_mass=dry_mass, charge=charge, ratios=ratios,
+                      additional_properties=additional_properties,
+                      )
 
     def _init_ws(self, pH=7., SAlk=None, COD=None, BOD=None,
                  uBOD=None, ThOD=None, cnBOD=None,
                  TC=None, TOC=None, TN=None, TKN=None,
                  TP=None, TK=None, TMg=None, TCa=None,
-                 dry_mass=None, charge=None, ratios=None):
+                 dry_mass=None, charge=None, ratios=None,
+                 additional_properties={}):
 
         self._pH = pH
         self._SAlk = SAlk
@@ -306,6 +317,7 @@ def _init_ws(self, pH=7., SAlk=None, COD=None, BOD=None,
         self._ratios = ratios
         self._state = None
         self._dstate = None
+        self.additional_properties = {}
 
     @staticmethod
     def from_stream(stream, ID='', **kwargs):
@@ -2241,6 +2253,11 @@ def dry_mass(self):
     def density(self):
         '''[float] Density of the stream, in mg/L (kg/m3).'''
         return 0.
+
+    @property
+    def additional_property(self):
+        '''[dict] Additional properties (e.g., turbidity).'''
+        return {}
 
     #!!! Keep this up-to-date with WasteStream
     # @property

qsdsan/data/process_data/_madm1.tsv

@@ -0,0 +1,47 @@
+	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_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
+hydrolysis_carbs	1									?	?	?		-1																															
+hydrolysis_proteins		1								?	?	?			-1																														
+hydrolysis_lipids	1-f_fa_li		f_fa_li							?	?	?				-1																													
+uptake_sugars	-1				(1-Y_su)*f_bu_su	(1-Y_su)*f_pro_su	(1-Y_su)*f_ac_su	(1-Y_su)*f_h2_su		?	?	?					Y_su																												
+uptake_amino_acids		-1		(1-Y_aa)*f_va_aa	(1-Y_aa)*f_bu_aa	(1-Y_aa)*f_pro_aa	(1-Y_aa)*f_ac_aa	(1-Y_aa)*f_h2_aa		?	?	?						Y_aa																											
+uptake_LCFA			-1				(1-Y_fa)*f_ac_fa	(1-Y_fa)*f_h2_fa		?	?	?							Y_fa																										
+uptake_valerate				-1		(1-Y_c4)*f_pro_va	(1-Y_c4)*f_ac_va	(1-Y_c4)*f_h2_va		?	?	?								Y_c4																									
+uptake_butyrate					-1		(1-Y_c4)*f_ac_bu	(1-Y_c4)*f_h2_bu		?	?	?								Y_c4																									
+uptake_propionate						-1	(1-Y_pro)*f_ac_pro	(1-Y_pro)*f_h2_pro		?	?	?									Y_pro																								
+uptake_acetate							-1		1-Y_ac	?	?	?										Y_ac																							
+uptake_h2								-1	1-Y_h2	?	?	?											Y_h2																						
+decay_Xsu										?	?	?	f_sI_xb	f_ch_xb	f_pr_xb	f_li_xb	-1							f_xI_xb																					
+decay_Xaa										?	?	?	f_sI_xb	f_ch_xb	f_pr_xb	f_li_xb		-1						f_xI_xb																					
+decay_Xfa										?	?	?	f_sI_xb	f_ch_xb	f_pr_xb	f_li_xb			-1					f_xI_xb																					
+decay_Xc4										?	?	?	f_sI_xb	f_ch_xb	f_pr_xb	f_li_xb				-1				f_xI_xb																					
+decay_Xpro										?	?	?	f_sI_xb	f_ch_xb	f_pr_xb	f_li_xb					-1			f_xI_xb																					
+decay_Xac										?	?	?	f_sI_xb	f_ch_xb	f_pr_xb	f_li_xb						-1		f_xI_xb																					
+decay_Xh2										?	?	?	f_sI_xb	f_ch_xb	f_pr_xb	f_li_xb							-1	f_xI_xb																					
+storage_PHA_valerate				-1						?	?	?													1	-Y_PO4		Y_PO4*K_XPP	Y_PO4*Mg_XPP																
+storage_PHA_butyrate					-1					?	?	?													1	-Y_PO4		Y_PO4*K_XPP	Y_PO4*Mg_XPP																
+storage_PHA_propionate						-1				?	?	?													1	-Y_PO4		Y_PO4*K_XPP	Y_PO4*Mg_XPP																
+storage_PHA_actate							-1			?	?	?													1	-Y_PO4		Y_PO4*K_XPP	Y_PO4*Mg_XPP																
+lysis_XPAO										?	?	?	f_sI_xb	f_ch_xb	f_pr_xb	f_li_xb								f_xI_xb			-1																		
+lysis_XPP										?	?	?														-1		K_XPP	Mg_XPP																
+lysis_XPHA				f_va_pha	f_bu_pha	f_pro_pha	f_ac_pha			?	?	?													-1																				
+growth_SRB_h2								-1		?	?	?																		-1*(1-Y_hSRB)*i_mass_IS*(MW_S0/MW_IS)	1-Y_hSRB	Y_hSRB													
+decay_XhSRB										?	?	?	f_sI_xb	f_ch_xb	f_pr_xb	f_li_xb								f_xI_xb								-1													
+growth_SRB_acetate							-1			?	?	?																		-1*(1-Y_aSRB)*i_mass_IS*(MW_S0/MW_IS)	1-Y_aSRB		Y_aSRB												
+decay_XaSRB										?	?	?	f_sI_xb	f_ch_xb	f_pr_xb	f_li_xb								f_xI_xb									-1												
+growth_SRB_propionate						-1	(1-Y_pSRB)*f_ac_pro			?	?	?																		-1*(1-Y_pSRB)*f_is_pro*i_mass_IS*(MW_S0/MW_IS)	(1-Y_pSRB)*f_is_pro			Y_pSRB											
+decay_XpSRB										?	?	?	f_sI_xb	f_ch_xb	f_pr_xb	f_li_xb								f_xI_xb										-1											
+growth_SRB_butyrate					-1		(1-Y_c4SRB)*f_ac_bu			?	?	?																		-1*(1-Y_c4SRB)*f_is_bu*i_mass_IS*(MW_S0/MW_IS)	(1-Y_c4SRB)*f_is_bu				Y_c4SRB										
+growth_SRB_valerate				-1		(1-Y_c4SRB)*f_pro_va	(1-Y_c4SRB)*f_ac_va			?	?	?																		-1*(1-Y_c4SRB)*f_is_va*i_mass_IS*(MW_S0/MW_IS)	(1-Y_c4SRB)*f_is_va				Y_c4SRB										
+decay_Xc4SRB										?	?	?	f_sI_xb	f_ch_xb	f_pr_xb	f_li_xb								f_xI_xb											-1										
+reduction_HFO_H_h2								-1																														1	-1*i_mass_Fe2						
+reduction_HFO_L_h2								-1																														1		-1*i_mass_Fe2					
+reduction_HFO_H_IS																															-1					(MW_S0/i_mass_S0)/(MW_IS/i_mass_IS)		1-(MW_S0/i_mass_S0)/(MW_IS/i_mass_IS)	i_mass_Fe2*((MW_S0/i_mass_S0)/(MW_IS/i_mass_IS)-1)						
+reduction_HFO_L_IS																															-1					(MW_S0/i_mass_S0)/(MW_IS/i_mass_IS)		1-(MW_S0/i_mass_S0)/(MW_IS/i_mass_IS)		i_mass_Fe2*((MW_S0/i_mass_S0)/(MW_IS/i_mass_IS)-1)					
+aging_HFO_H																																							-1	1					
+aging_HFO_L																																								-1	1				
+fast_P_binding												?																											-1			1			
+slow_P_sorption												?																												-1			1		
+aging_HFO_HP																																										-1		1	
+aging_HFO_LP																																											-1		1
+dissolution_HFO_HP												?																											1			-1			
+dissolution_HFO_LP												?																												1			-1		

qsdsan/processes/__init__.py

@@ -17,6 +17,7 @@
 from ._asm1 import *
 from ._asm2d import *
 from ._adm1 import *
+from ._madm1 import *
 from ._decay import *
 from ._kinetic_reaction import *
 from ._pm2 import *
@@ -26,6 +27,7 @@
     _asm1,
     _asm2d,
     _adm1,
+    _madm1,
     _decay,
     _kinetic_reaction,
     )
@@ -35,6 +37,7 @@
     *_asm1.__all__,
     *_asm2d.__all__,
     *_adm1.__all__,
+    *_madm1.__all__,
     *_decay.__all__,
     *_kinetic_reaction.__all__,
     *_pm2.__all__,

qsdsan/processes/_adm1.py

@@ -190,6 +190,9 @@ def mass2mol_conversion(cmps):
 #     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):
@@ -267,10 +270,24 @@ def rhos_adm1(state_arr, params):
     weak_acids = cmps_in_M[[24, 25, 10, 9, 6, 5, 4, 3]]
 
     T_op = state_arr[-1]
+    if T_op == T_base:
+        Ka = Kab
+        KH = KHb / unit_conversion[7:10]
+    else:
+        T_temp = params.pop('T_op', None)
+        if T_op == T_temp:
+            params['T_op'] = T_op
+            Ka = params['Ka']
+            KH = params['KH']
+        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:10]
+
     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]
+    # 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)
@@ -279,12 +296,12 @@ def rhos_adm1(state_arr, params):
     if S_bu > 0: rhos[8] *= 1/(1+S_va/S_bu)
 
     h = brenth(acid_base_rxn, 1e-14, 1.0,
-            args=(weak_acids, Kas),
+            args=(weak_acids, Ka),
             xtol=1e-12, maxiter=100)
     # h = 10**(-7.46)
 
-    nh3 = Kas[1] * weak_acids[2] / (Kas[1] + h)
-    co2 = weak_acids[3] - Kas[2] * weak_acids[3] / (Kas[2] + h)
+    nh3 = Ka[1] * weak_acids[2] / (Ka[1] + h)
+    co2 = weak_acids[3] - Ka[2] * weak_acids[3] / (Ka[2] + h)
     biogas_S[-1] = co2 / unit_conversion[9]
 
     Iph = Hill_inhibit(h, pH_ULs, pH_LLs)