Skip to content

Commit

Permalink
improvements to phenomena-oriented simulation and add tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@yoelcortes
yoelcortes committed Jan 30, 2024
1 parent eed7cfc commit 4c65749
Show file tree
Hide file tree
Showing 5 changed files with 139 additions and 128 deletions.
3 changes: 1 addition & 2 deletions biosteam/_network.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -256,7 +256,6 @@ def sort(self, ends):
path_sources = [PathSource(i, ends) for i in self.path]
N = len(path_sources)
if not N: return
has_recycle = bool(self.recycle)
for _ in range(N * N):
stop = True
for i in range(N - 1):
Expand All @@ -281,7 +280,7 @@ def sort(self, ends):
else:
path_sources.remove(downstream)
path_sources.insert(i, downstream)
upstream, downstream = downstream, upstream
upstream = downstream
stop = False
break
if stop: break
Expand Down
20 changes: 13 additions & 7 deletions biosteam/_system.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@ def __init__(self, variable, units=None):
def append(self, unit):
A = self.A; b = self.b
self.units.append(unit)
zero = np.array(0)
# zero = np.array(0)
for coefficients, value in unit._create_linear_equations(self.variable):
# if all([(i == zero).all() for i in coefficients.values()]):
# continue
Expand Down Expand Up @@ -120,11 +120,16 @@ def solve(self):
values.append(value)
return objs, np.array(values)
A, objs = dictionaries2array(self.A)
# (A == 0).all(axis=1)
# rows = A.any(axis=0)
# cols = A.any(axis=1)
# b = [j for (i, j) in zip(rows, b)]
b = np.array(b).T
# A = A[rows, ...][:, cols, ...]
# objs = [j for (i, j) in zip(cols, objs) if i]
try:
values = solve(A, b).T
except Exception as e:
breakpoint()
for i in self.A:
print('--')
for i, j in i.items():
Expand Down Expand Up @@ -677,7 +682,7 @@ class System:
available_methods: Methods[str, tuple[Callable, bool, dict]] = Methods()

#: Variable solution priority for phenomena oriented simulation.
variable_priority: list[str] = [('mol', 'K-pseudo'), 'K', 'T', 'L', 'B']
variable_priority: list[str] = ['mol', ('mol-LLE', 'K-pseudo'), 'K', 'T', 'L', 'B']

@classmethod
def register_method(cls, name, solver, conditional=False, **kwargs):
Expand Down Expand Up @@ -2207,14 +2212,15 @@ def run(self):
if algorithm == 'Sequential modular':
self.run_sequential_modular()
elif algorithm == 'Phenomena oriented':
if self._iter == 0:
self.run_sequential_modular()
if all([i.isempty() for i in self.get_all_recycles()]):
for i in self.unit_path: i.run()
else:
try:
self.run_phenomena()
except:
print('FAILED PHENOMENA-ORIENTED SIMULATION')
self.run_sequential_modular()
warn('phenomena-oriented simulation failed; '
'attempting one sequential-modular loop', RuntimeWarning)
for i in self.unit_path: i.run()
else:
raise RuntimeError(f'unknown algorithm {algorithm!r}')

Expand Down
14 changes: 12 additions & 2 deletions biosteam/units/phase_equilibrium.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -530,13 +530,23 @@ def _run_decoupled_Kgamma(self, P=None): # Psuedo-equilibrium
f_gamma, IDs, index = self._get_activity_model()
T = self.T
x = bottom.mol[index]
x /= x.sum()
x_sum = x.sum()
if x_sum:
x /= x_sum
else:
x = np.ones(x.size) / x.size
gamma_x = f_gamma(x, T)
gamma_y = self.gamma_y
if gamma_y is None or gamma_y.size != len(index):
try:
init_gamma = gamma_y is None or gamma_y.size != len(index)
except:
init_gamma = True
if init_gamma:
y = top.mol[index]
y /= y.sum()
self.gamma_y = gamma_y = f_gamma(y, T)


K = self.K
K = gamma_x / gamma_y
y = K * x
Expand Down
228 changes: 112 additions & 116 deletions tests/test_phenomena_oriented_simulation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,7 +104,7 @@ def system(ins, outs):
relative_molar_tolerance=1e-6,
method='wegstein')

for i in range(1):
for i in range(2):
dp_sys.simulate()
dp_sys.run_phenomena()
for i in range(2):
Expand All @@ -117,6 +117,9 @@ def system(ins, outs):
assert_allclose(actual, value, rtol=0.001)

def test_simple_acetic_acid_separation_with_recycle():
import biosteam as bst
import numpy as np
from numpy.testing import assert_allclose
bst.settings.set_thermo(['Water', 'AceticAcid', 'EthylAcetate'], cache=True)
solvent_feed_ratio = 1
@bst.SystemFactory
Expand Down Expand Up @@ -194,130 +197,123 @@ def fresh_solvent_flow_rate():
for i in range(1): sm_sys.simulate()
t1 = time.toc()

assert t0 < 0.8 * t1, 'phenomena-oriented simulation speed should be faster than sequential-modular'
assert t0 < 0.75 * t1, 'phenomena-oriented simulation speed should be faster than sequential-modular'

for s_sm, s_dp in zip(sm_sys.streams, dp_sys.streams):
actual = s_sm.mol
value = s_dp.mol
assert_allclose(actual, value, rtol=0.01)

# def test_complex_acetic_acid_separation_system():
# import biosteam as bst
# import numpy as np
# bst.settings.set_thermo(['Water', 'AceticAcid', 'EthylAcetate'], cache=True)
# def create_system(alg):
# solvent_feed_ratio = 1
# chemicals = bst.settings.chemicals
# acetic_acid_broth = bst.Stream(
# ID='acetic_acid_broth', AceticAcid=6660, Water=43600, units='kg/hr'
# )
# ethyl_acetate = bst.Stream(
# ID='fresh_solvent', EthylAcetate=15000, units='kg/hr'
# )
# glacial_acetic_acid = bst.Stream('glacial_acetic_acid')
# wastewater = bst.Stream('wastewater')
# solvent_recycle = bst.Stream('solvent_rich')
# with bst.System(algorithm=alg) as sys:
# @ethyl_acetate.equation('mol')
# def fresh_solvent_flow_rate():
# f = np.ones(chemicals.size)
# r = np.zeros(chemicals.size)
# v = r.copy()
# index = chemicals.index('EthylAcetate')
# r[index] = 1
# v[index] = solvent_feed_ratio * acetic_acid_broth.F_mass / chemicals.EthylAcetate.MW
# return (
# {ethyl_acetate: f,
# solvent_recycle: r},
# v
# )
# extractor = bst.MultiStageMixerSettlers(
# 'extractor',
# ins=(acetic_acid_broth, ethyl_acetate, solvent_recycle),
# outs=('extract', 'raffinate'),
# top_chemical='EthylAcetate',
# feed_stages=(0, -1, -1),
# N_stages=5,
# collapsed_init=False,
# use_cache=True,
# )
def test_complex_acetic_acid_separation_system():
import biosteam as bst
import numpy as np
bst.settings.set_thermo(['Water', 'AceticAcid', 'EthylAcetate'], cache=True)
def create_system(alg):
solvent_feed_ratio = 1
chemicals = bst.settings.chemicals
acetic_acid_broth = bst.Stream(
ID='acetic_acid_broth', AceticAcid=6660, Water=43600, units='kg/hr'
)
ethyl_acetate = bst.Stream(
ID='fresh_solvent', EthylAcetate=15000, units='kg/hr'
)
glacial_acetic_acid = bst.Stream('glacial_acetic_acid')
wastewater = bst.Stream('wastewater')
solvent_recycle = bst.Stream('solvent_rich')
with bst.System(algorithm=alg) as sys:
@ethyl_acetate.equation('mol')
def fresh_solvent_flow_rate():
f = np.ones(chemicals.size)
r = np.zeros(chemicals.size)
v = r.copy()
index = chemicals.index('EthylAcetate')
r[index] = 1
v[index] = solvent_feed_ratio * acetic_acid_broth.F_mass / chemicals.EthylAcetate.MW
return (
{ethyl_acetate: f,
solvent_recycle: r},
v
)
extractor = bst.MultiStageMixerSettlers(
'extractor',
ins=(acetic_acid_broth, ethyl_acetate, solvent_recycle),
outs=('extract', 'raffinate'),
top_chemical='EthylAcetate',
feed_stages=(0, -1, -1),
N_stages=5,
collapsed_init=False,
use_cache=True,
)

# @extractor.add_specification(run=True)
# def adjust_fresh_solvent_flow_rate():
# broth = acetic_acid_broth.F_mass
# EtAc_recycle = solvent_recycle.imass['EthylAcetate']
# ethyl_acetate.imass['EthylAcetate'] = max(
# 0, broth * solvent_feed_ratio - EtAc_recycle
# )
# reflux = bst.Stream()
# ED = bst.MESHDistillation(
# 'extract_distiller',
# ins=[extractor.extract, reflux],
# outs=['vapor', glacial_acetic_acid],
# LHK=('EthylAcetate', 'AceticAcid'),
# N_stages=8,
# feed_stages=(4, 0),
# reflux=None,
# boilup=2.,
# )
# condenser = bst.StageEquilibrium(
# 'condenser',
# ins=ED-0,
# outs=('', 'condensate'),
# phases=('g', 'l'),
# B=0,
# )
# water_rich = bst.Stream('water_rich')
# distillate = bst.Stream('distillate')
# settler = bst.StageEquilibrium(
# 'settler',
# ins=(condenser-1, distillate),
# outs=(solvent_recycle, water_rich, reflux),
# phases=('L', 'l'),
# top_chemical='EthylAcetate',
# top_split=0.5,
# )
# RD = bst.MESHDistillation(
# 'raffinate_distiller',
# LHK=('EthylAcetate', 'Water'),
# ins=[extractor.raffinate, water_rich],
# outs=['', wastewater, distillate],
# full_condenser=True,
# N_stages=8,
# feed_stages=(3, 3),
# reflux=1.5,
# boilup=2.,
# )
# return sys
# time = bst.TicToc()
# sys_init = create_system('sequential modular')
# time.tic()
# sys_init.set_tolerance(rmol=0.1, mol=1e-6, subsystems=True)
# sys_init.simulate()
# t_init = time.toc()

# time.tic()
# po = create_system('sequential modular')
# po.set_tolerance(rmol=0.1, mol=1e-6, subsystems=True)
# po.simulate()
# po.set_tolerance(rmol=1e-3, mol=1e-6, subsystems=True)
# po.algorithm = 'phenomena oriented'
# po.method = 'fixed-point'
# po.simulate()
# t_phenomena = time.toc()
@extractor.add_specification(run=True)
def adjust_fresh_solvent_flow_rate():
broth = acetic_acid_broth.F_mass
EtAc_recycle = solvent_recycle.imass['EthylAcetate']
ethyl_acetate.imass['EthylAcetate'] = max(
0, broth * solvent_feed_ratio - EtAc_recycle
)
reflux = bst.Stream()
ED = bst.MESHDistillation(
'extract_distiller',
ins=[extractor.extract, reflux],
outs=['vapor', glacial_acetic_acid],
LHK=('EthylAcetate', 'AceticAcid'),
N_stages=8,
feed_stages=(4, 0),
reflux=None,
boilup=2.,
)
condenser = bst.StageEquilibrium(
'condenser',
ins=ED-0,
outs=('', 'condensate'),
phases=('g', 'l'),
B=0,
)
water_rich = bst.Stream('water_rich')
distillate = bst.Stream('distillate')
settler = bst.StageEquilibrium(
'settler',
ins=(condenser-1, distillate),
outs=(solvent_recycle, water_rich, reflux),
phases=('L', 'l'),
top_chemical='EthylAcetate',
top_split=0.5,
)
RD = bst.MESHDistillation(
'raffinate_distiller',
LHK=('EthylAcetate', 'Water'),
ins=[extractor.raffinate, water_rich],
outs=['', wastewater, distillate],
full_condenser=True,
N_stages=8,
feed_stages=(3, 3),
reflux=1.5,
boilup=2.,
)
return sys
time = bst.TicToc()

time.tic()
sm = create_system('sequential modular')
sm.set_tolerance(rmol=1e-6, mol=1e-6, subsystems=True, method='fixed-point')
sm.simulate()
t_sequential = time.toc()

# time.tic()
# sm = create_system('sequential modular')
# sm.set_tolerance(rmol=1e-3, mol=1e-6, subsystems=True)
# sm.simulate()
# t_sequential = time.toc()
time.tic()
po = create_system('phenomena oriented')
po.set_tolerance(rmol=1e-6, mol=1e-6,
subsystems=True,
method='fixed-point')
po.simulate()
t_phenomena = time.toc()

# print(t_init, t_phenomena, t_sequential)
assert t_phenomena < 0.2 * t_sequential

# for s_sm, s_dp in zip(sm.streams, po.streams):
# actual = s_sm.mol
# value = s_dp.mol
# assert_allclose(actual, value, rtol=0.01)
for s_sm, s_dp in zip(sm.streams, po.streams):
actual = s_sm.mol
value = s_dp.mol
assert_allclose(actual, value, rtol=0.01)

# def test_complex_acetic_acid_separation_system():
# import biosteam as bst
Expand Down Expand Up @@ -460,4 +456,4 @@ def fresh_solvent_flow_rate():
test_trivial_distillation_case()
test_simple_acetic_acid_separation_no_recycle()
test_simple_acetic_acid_separation_with_recycle()
# test_complex_acetic_acid_separation_system()
test_complex_acetic_acid_separation_system()
2 changes: 1 addition & 1 deletion thermosteam
Submodule thermosteam updated 1 files
+2 −2 thermosteam/equilibrium/activity_coefficients.py

0 comments on commit 4c65749

Please sign in to comment.