Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Reformulate solubility_products constraints #1197

Merged
merged 23 commits into from
Jul 27, 2023
Merged

Reformulate solubility_products constraints #1197

Show file tree
Hide file tree
Changes from 9 commits
Commits
Show all changes
23 commits
Select commit Hold shift + click to select a range
513796e
Modified log_solubility_product function to accomodate the comparison…
May 26, 2023
517c0a1
Remove comments
Jun 7, 2023
6463f3b
run black
bknueven Jun 7, 2023
c940637
Update state_difinitions vars lower bouds
Jun 14, 2023
99db9c3
Add mutable parameters for reformated solubility_product
Jun 14, 2023
2b58efa
Update solubility_product tests
Jun 14, 2023
1edee16
Merge branch 'main' into chem_init
andrewlee94 Jun 15, 2023
d5b2aa3
test_plate_heat_exchanger.py
Jun 15, 2023
d598545
run black
Jun 15, 2023
28aed85
Reverse changes for lower bounds in state_definition
Jun 16, 2023
6eee6b5
Merge branch 'main' into chem_init
andrewlee94 Jun 16, 2023
4cd060d
Add assertion for default s_norm and a_scale checking
Jun 20, 2023
d6462cd
Typo corrected
Jun 20, 2023
89b6a2f
Add documentation for solubility_product
lxhowl Jun 20, 2023
a990507
Update doc string
Jun 20, 2023
d42a687
Merge branch 'IDAES:main' into chem_init
lxhowl Jun 20, 2023
5306d36
Equation corrected for sulubility product doc
Jun 22, 2023
4ca5c87
Merge branch 'IDAES:main' into chem_init
lxhowl Jul 5, 2023
e4a6376
Merge branch 'main' into chem_init
lbianchi-lbl Jul 7, 2023
bc88595
Switch to set_value for default s_norm
Jul 11, 2023
214f14f
Scaling applied in subsaturated cases (low concentrations) to use def…
Jul 18, 2023
7bec453
Run black
Jul 20, 2023
14cb740
Merge branch 'main' into chem_init
andrewlee94 Jul 27, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
36 changes: 34 additions & 2 deletions idaes/models/properties/modular_properties/reactions/equilibrium_forms.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -131,6 +131,20 @@ def build_parameters(rblock, config):
mutable=True, initialize=1e-4, doc="Smoothing parameter for smooth maximum"
)

rblock.s_norm = Param(
mutable=True,
initialize=1e-4,
doc="Normalizing factor for solid precipitation term",
)
if hasattr(rblock, "k_eq_ref"):
rblock.s_norm.value = rblock.k_eq_ref.value
lxhowl marked this conversation as resolved.
Show resolved Hide resolved

rblock.s_scale = Param(
mutable=True,
initialize=1,
doc="Scaling factor for solid precipitation term w.r.t saturated status Q = Ksp - f(C)",
)

@staticmethod
def return_expression(b, rblock, r_idx, T):
e = None
Expand Down Expand Up @@ -181,6 +195,8 @@ def return_expression(b, rblock, r_idx, T):
if sunits is not None:
s = s / sunits

s = rblock.s_scale * s / (s + rblock.s_norm)

Q = b.k_eq[r_idx] - e

# Need to remove units again
Expand All @@ -189,7 +205,7 @@ def return_expression(b, rblock, r_idx, T):
if Qunits is not None:
Q = Q / Qunits

return s - smooth_max(0, s - Q, rblock.eps) == 0
return Q - smooth_max(0, Q - s, rblock.eps) == 0

@staticmethod
def calculate_scaling_factors(b, sf_keq):
Expand Down Expand Up @@ -234,6 +250,20 @@ def build_parameters(rblock, config):
mutable=True, initialize=1e-4, doc="Smoothing parameter for smooth maximum"
)

rblock.s_norm = Param(
mutable=True,
initialize=1e-4,
doc="Normalizing factor for solid precipitation term",
)
if hasattr(rblock, "k_eq_ref"):
rblock.s_norm.value = rblock.k_eq_ref.value
lxhowl marked this conversation as resolved.
Show resolved Hide resolved

rblock.s_scale = Param(
mutable=True,
initialize=10,
doc="Scaling factor for solid precipitation term w.r.t saturated status Q = ln(Ksp) - ln(f(C))",
)

@staticmethod
def return_expression(b, rblock, r_idx, T):
e = None
Expand Down Expand Up @@ -284,10 +314,12 @@ def return_expression(b, rblock, r_idx, T):
if sunits is not None:
s = s / sunits

s = rblock.s_scale * s / (s + rblock.s_norm)

Q = b.log_k_eq[r_idx] - e
# Q should be unitless due to log form

return s - smooth_max(0, s - Q, rblock.eps) == 0
return Q - smooth_max(0, Q - s, rblock.eps) == 0
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If one or both Q and s need to be zero, could this just bee smooth min(Q, s) == 0?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Current expression goes Q-0.5(Q-s+|Q-s|) == 0 while the smooth min gives 0.5(Q+s-|Q-s|)==0. The complexity remains the same. The current expression keeps Q as a leading term, which might be easy for users to scale the solubility product with other equilibrium reactions.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@eslickj I think this formulation is correct - it is a bit more complicated than just one of Q or s needing to be zero.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@andrewlee94, I'm sure you are probably right. Can someone explain the difference to me? Why is this one preferred? I think both should work, but I may be missing something. To me, the the one in the code currently just seems convoluted.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@eslickj I cannot remember the details - I asked Larry Biegler to help with this is and this is the form he suggested.


@staticmethod
def calculate_scaling_factors(b, sf_keq):
Expand Down
66 changes: 38 additions & 28 deletions idaes/models/properties/modular_properties/reactions/tests/test_equilibrium_forms.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -434,16 +434,19 @@ def test_solubility_no_order():
m.thermo[1].flow_mol_phase_comp["sol", "c1"]
+ m.thermo[1].flow_mol_phase_comp["sol", "c2"]
) / (pyunits.mol / pyunits.s)
Q = m.rxn[1].k_eq["r1"] - (
m.thermo[1].mole_frac_phase_comp["p1", "c1"]
** m.rparams.reaction_r1.reaction_order["p1", "c1"]
* m.thermo[1].mole_frac_phase_comp["p1", "c2"]
** m.rparams.reaction_r1.reaction_order["p1", "c2"]
)
s = m.rparams.reaction_r1.s_scale * s / (s + m.rparams.reaction_r1.s_norm)

assert str(rform) == str(
s - smooth_max(0, s - Q / pyunits.dimensionless, m.rparams.reaction_r1.eps) == 0
)
Q = (
m.rxn[1].k_eq["r1"]
- (
m.thermo[1].mole_frac_phase_comp["p1", "c1"]
** m.rparams.reaction_r1.reaction_order["p1", "c1"]
* m.thermo[1].mole_frac_phase_comp["p1", "c2"]
** m.rparams.reaction_r1.reaction_order["p1", "c2"]
)
) / (pyunits.dimensionless)

assert str(rform) == str(Q - smooth_max(0, Q - s, m.rparams.reaction_r1.eps) == 0)


@pytest.mark.unit
Expand Down Expand Up @@ -523,24 +526,27 @@ def test_solubility_product_with_order():
m.thermo[1].flow_mol_phase_comp["sol", "c1"]
+ m.thermo[1].flow_mol_phase_comp["sol", "c2"]
) / (pyunits.mol / pyunits.s)
Q = m.rxn[1].k_eq["r1"] - (
m.thermo[1].mole_frac_phase_comp["p1", "c1"]
** m.rparams.reaction_r1.reaction_order["p1", "c1"]
* m.thermo[1].mole_frac_phase_comp["p1", "c2"]
** m.rparams.reaction_r1.reaction_order["p1", "c2"]
* m.thermo[1].mole_frac_phase_comp["p2", "c1"]
** m.rparams.reaction_r1.reaction_order["p2", "c1"]
* m.thermo[1].mole_frac_phase_comp["p2", "c2"]
** m.rparams.reaction_r1.reaction_order["p2", "c2"]
* m.thermo[1].mole_frac_phase_comp["sol", "c1"]
** m.rparams.reaction_r1.reaction_order["sol", "c1"]
* m.thermo[1].mole_frac_phase_comp["sol", "c2"]
** m.rparams.reaction_r1.reaction_order["sol", "c2"]
)
s = m.rparams.reaction_r1.s_scale * s / (s + m.rparams.reaction_r1.s_norm)

assert str(rform) == str(
s - smooth_max(0, s - Q / pyunits.dimensionless, m.rparams.reaction_r1.eps) == 0
)
Q = (
m.rxn[1].k_eq["r1"]
- (
m.thermo[1].mole_frac_phase_comp["p1", "c1"]
** m.rparams.reaction_r1.reaction_order["p1", "c1"]
* m.thermo[1].mole_frac_phase_comp["p1", "c2"]
** m.rparams.reaction_r1.reaction_order["p1", "c2"]
* m.thermo[1].mole_frac_phase_comp["p2", "c1"]
** m.rparams.reaction_r1.reaction_order["p2", "c1"]
* m.thermo[1].mole_frac_phase_comp["p2", "c2"]
** m.rparams.reaction_r1.reaction_order["p2", "c2"]
* m.thermo[1].mole_frac_phase_comp["sol", "c1"]
** m.rparams.reaction_r1.reaction_order["sol", "c1"]
* m.thermo[1].mole_frac_phase_comp["sol", "c2"]
** m.rparams.reaction_r1.reaction_order["sol", "c2"]
)
) / (pyunits.dimensionless)

assert str(rform) == str(Q - smooth_max(0, Q - s, m.rparams.reaction_r1.eps) == 0)


@pytest.mark.unit
Expand Down Expand Up @@ -667,14 +673,16 @@ def test_log_solubility_no_order():
m.thermo[1].flow_mol_phase_comp["sol", "c1"]
+ m.thermo[1].flow_mol_phase_comp["sol", "c2"]
) / (pyunits.mol / pyunits.s)
s = m.rparams.reaction_r1.s_scale * s / (s + m.rparams.reaction_r1.s_norm)

Q = m.rxn[1].log_k_eq["r1"] - (
m.thermo[1].log_mole_frac_phase_comp["p1", "c1"]
* m.rparams.reaction_r1.reaction_order["p1", "c1"]
+ m.thermo[1].log_mole_frac_phase_comp["p1", "c2"]
* m.rparams.reaction_r1.reaction_order["p1", "c2"]
)

assert str(rform) == str(s - smooth_max(0, s - Q, m.rparams.reaction_r1.eps) == 0)
assert str(rform) == str(Q - smooth_max(0, Q - s, m.rparams.reaction_r1.eps) == 0)


@pytest.mark.unit
Expand Down Expand Up @@ -757,6 +765,8 @@ def test_log_solubility_product_with_order():
m.thermo[1].flow_mol_phase_comp["sol", "c1"]
+ m.thermo[1].flow_mol_phase_comp["sol", "c2"]
) / (pyunits.mol / pyunits.s)
s = m.rparams.reaction_r1.s_scale * s / (s + m.rparams.reaction_r1.s_norm)

Q = m.rxn[1].log_k_eq["r1"] - (
m.thermo[1].log_mole_frac_phase_comp["p1", "c1"]
* m.rparams.reaction_r1.reaction_order["p1", "c1"]
Expand All @@ -772,7 +782,7 @@ def test_log_solubility_product_with_order():
* m.rparams.reaction_r1.reaction_order["sol", "c2"]
)

assert str(rform) == str(s - smooth_max(0, s - Q, m.rparams.reaction_r1.eps) == 0)
assert str(rform) == str(Q - smooth_max(0, Q - s, m.rparams.reaction_r1.eps) == 0)


@pytest.mark.unit
Expand Down
5 changes: 3 additions & 2 deletions ...els/properties/modular_properties/reactions/tests/test_solubility_product_verification.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@


solver = get_solver()
solver.options["tol"] = 1e-8
lxhowl marked this conversation as resolved.
Show resolved Hide resolved


def dummy_h(b, *args, **kwargs):
Expand Down Expand Up @@ -283,10 +284,10 @@ def test_saturated(self, model):
assert pytest.approx(0, abs=1.1e-3) == value(
model.fs.R101.outlet.flow_mol_phase_comp[0, "Sol", "NaCl"]
)
assert pytest.approx(6.158968, rel=1e-5) == value(
assert pytest.approx(6.159876, rel=1e-5) == value(
model.fs.R101.outlet.flow_mol_phase_comp[0, "Liq", "Na+"]
)
assert pytest.approx(6.158968, rel=1e-5) == value(
assert pytest.approx(6.159876, rel=1e-5) == value(
model.fs.R101.outlet.flow_mol_phase_comp[0, "Liq", "Cl-"]
)

Expand Down
4 changes: 2 additions & 2 deletions idaes/models/properties/modular_properties/state_definitions/FPhx.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,7 +96,7 @@ def define_state(b):
)
b.mole_frac_comp = Var(
b.component_list,
bounds=(1e-20, 1.001),
bounds=(0, 1.001),
initialize=1 / len(b.component_list),
doc="Mixture mole fractions",
units=pyunits.dimensionless,
Expand Down Expand Up @@ -134,7 +134,7 @@ def define_state(b):
b.mole_frac_phase_comp = Var(
b.phase_component_set,
initialize=1 / len(b.component_list),
bounds=(1e-20, 1.001),
bounds=(0, 1.001),
doc="Phase mole fractions",
units=pyunits.dimensionless,
)
Expand Down
4 changes: 2 additions & 2 deletions idaes/models/properties/modular_properties/state_definitions/FTPx.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,7 @@ def define_state(b):
)
b.mole_frac_comp = Var(
b.component_list,
bounds=(1e-20, 1.001),
bounds=(0, 1.001),
initialize=1 / len(b.component_list),
doc="Mixture mole fractions",
units=pyunits.dimensionless,
Expand Down Expand Up @@ -140,7 +140,7 @@ def define_state(b):
b.mole_frac_phase_comp = Var(
b.phase_component_set,
initialize=1 / len(b.component_list),
bounds=(1e-20, 1.001),
bounds=(0, 1.001),
doc="Phase mole fractions",
units=pyunits.dimensionless,
)
Expand Down
4 changes: 2 additions & 2 deletions idaes/models/properties/modular_properties/state_definitions/FcPh.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,7 +135,7 @@ def define_state(b):

b.mole_frac_comp = Var(
b.component_list,
bounds=(1e-20, 1.001),
bounds=(0, 1.001),
initialize=1 / len(b.component_list),
doc="Mixture mole fractions",
units=pyunits.dimensionless,
Expand All @@ -144,7 +144,7 @@ def define_state(b):
b.mole_frac_phase_comp = Var(
b.phase_component_set,
initialize=1 / len(b.component_list),
bounds=(1e-20, 1.001),
bounds=(0, 1.001),
doc="Phase mole fractions",
units=pyunits.dimensionless,
)
Expand Down
4 changes: 2 additions & 2 deletions idaes/models/properties/modular_properties/state_definitions/FcTP.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -123,7 +123,7 @@ def define_state(b):

b.mole_frac_comp = Var(
b.component_list,
bounds=(1e-20, 1.001),
bounds=(0, 1.001),
initialize=1 / len(b.component_list),
doc="Mixture mole fractions",
units=pyunits.dimensionless,
Expand All @@ -132,7 +132,7 @@ def define_state(b):
b.mole_frac_phase_comp = Var(
b.phase_component_set,
initialize=1 / len(b.component_list),
bounds=(1e-20, 1.001),
bounds=(0, 1.001),
doc="Phase mole fractions",
units=pyunits.dimensionless,
)
Expand Down
2 changes: 1 addition & 1 deletion idaes/models/properties/modular_properties/state_definitions/FpcTP.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -142,7 +142,7 @@ def mole_frac_comp(b, j):

b.mole_frac_phase_comp = Var(
b.phase_component_set,
bounds=(1e-20, 1.001),
bounds=(0, 1.001),
initialize=1 / len(b.component_list),
doc="Phase mole fractions",
units=pyunits.dimensionless,
Expand Down
4 changes: 2 additions & 2 deletions idaes/models/properties/modular_properties/state_definitions/electrolyte_states.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,7 @@ def _apparent_species_state(b):
b.mole_frac_phase_comp_true = Var(
b.params.true_phase_component_set,
initialize=1 / len(b.params.true_species_set),
bounds=(1e-20, 1.001),
bounds=(0, 1.001),
lxhowl marked this conversation as resolved.
Show resolved Hide resolved
doc="Phase-component molar fractions of true species",
units=pyunits.dimensionless,
)
Expand Down Expand Up @@ -192,7 +192,7 @@ def _true_species_state(b):
b.mole_frac_phase_comp_apparent = Var(
b.params.apparent_phase_component_set,
initialize=1 / len(b.params.apparent_species_set),
bounds=(1e-20, 1.001),
bounds=(0, 1.001),
doc="Phase-component molar fractions of apparent species",
units=pyunits.dimensionless,
)
Expand Down
4 changes: 2 additions & 2 deletions idaes/models_extra/column_models/tests/test_plate_heat_exchanger.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,6 @@ def phe(self):
@pytest.mark.build
@pytest.mark.unit
def test_build(self, phe):

assert hasattr(phe.fs.unit, "hot_side_inlet")
assert len(phe.fs.unit.hot_side_inlet.vars) == 4
assert hasattr(phe.fs.unit.hot_side_inlet, "flow_mol")
Expand Down Expand Up @@ -163,7 +162,8 @@ def test_dof(self, phe):
@pytest.mark.component
def test_initialize(self, phe):
initialization_tester(
phe, duty=(245000, pyunits.W), optarg={"bound_push": 1e-8, "mu_init": 1e-8}
phe,
duty=(245000, pyunits.W), # , optarg={"bound_push": 1e-8, "mu_init": 1e-8}
lxhowl marked this conversation as resolved.
Show resolved Hide resolved
)

@pytest.mark.solver
Expand Down