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.

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

Mea column improvements #1077

Draft
wants to merge 21 commits into
base: main
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
21 commits
Select commit Hold shift + click to select a range
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
13 changes: 12 additions & 1 deletion idaes/core/base/property_base.py
Original file line number Diff line number Diff line change
Expand Up @@ -761,11 +761,22 @@ def calculate_scaling_factors(self):
super().calculate_scaling_factors()
# Get scaling factor defaults, if no scaling factor set
for v in self.component_data_objects(
(Constraint, Var, Expression), descend_into=False
(Var, Expression), descend_into=False
):
if iscale.get_scaling_factor(v) is None: # don't replace if set
name = v.getname().split("[")[0]
index = v.index()
sf = self.config.parameters.get_default_scaling(name, index)
if sf is not None:
iscale.set_scaling_factor(v, sf)

for con in self.component_data_objects(Constraint, descend_into=False):
# if con.name == 'stripper_section.stripper.liquid_phase.properties[0.0,0.8].appr_to_true_species[Liq,CO2]':
Copy link
Member

Choose a reason for hiding this comment

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

This comment needs to be removed.

# import pdb; pdb.set_trace()
if iscale.get_constraint_transform_applied_scaling_factor(con) is None:
name = con.getname().split("[")[0]
index = con.index()
sf = self.config.parameters.get_default_scaling(name, index)
if sf is not None:
iscale.constraint_scaling_transform(con, sf, overwrite=False)
Copy link
Member

Choose a reason for hiding this comment

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

This is a fairly significant change to what this method does - we probably need to discuss the implications of this.

Copy link
Member

Choose a reason for hiding this comment

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

Thinking a bit more (and repeating a comment from Slack for posterity) - I think this is the wrong way to go.

IPOPT should be considered a special case, not the general case, in which case we need to set the scaling suffix for the constraints and then let the user decide what to do with them. IPOPT users can use the constraint_scaling_transformation function if they wish, but the generic workflow would use the Pyomo scaling transformation (which relies on having the suffixes).

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Removing the default status of IPOPT and constraint_scaling_transform in the long run is fine, but it's also a complete revolution in how scaling is treated in IDAES models, which make heavy use of constraint_scaling_transform. If this change is reverted, I'll just as likely use constraint_scaling_transform somewhere else to do the desired scaling: probably in electrolyte_states.py.

Copy link
Member

Choose a reason for hiding this comment

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

It is fine to still use constraint_scaling_transform, but it should not be done in the base classes (i.e. anything in idaes.core.base). Ideally, we would only do this in the model level calculate_scaling_factors method.

The point being this is not something we want baked in deep as it will be harder to change later - ideally we want to limit the use of this to the calculate_scaling_factors methods so that when we make the transition, we only need to change those (or if we are lucky, just switch to new methods and delete those as they are no longer needed).


Original file line number Diff line number Diff line change
Expand Up @@ -4151,7 +4151,7 @@ def _log_mole_frac_phase_comp(self):
try:
self.log_mole_frac_phase_comp = Var(
self.phase_component_set,
initialize=0,
initialize=-1,
bounds=(-100, log(1.001)),
units=pyunits.dimensionless,
doc="Log of mole fractions of component by phase",
Expand All @@ -4176,8 +4176,8 @@ def _log_mole_frac_phase_comp_apparent(self):
try:
self.log_mole_frac_phase_comp_apparent = Var(
self.params.apparent_phase_component_set,
initialize=0,
bounds=(-50, 0),
initialize=-1,
bounds=(-100, log(1.001)),
units=pyunits.dimensionless,
doc="Log of mole fractions of component by phase",
)
Expand All @@ -4201,8 +4201,8 @@ def _log_mole_frac_phase_comp_true(self):
try:
self.log_mole_frac_phase_comp_true = Var(
self.params.true_phase_component_set,
initialize=0,
bounds=(-50, 0),
initialize=-1,
bounds=(-100, log(1.001)),
units=pyunits.dimensionless,
doc="Log of mole fractions of component by phase",
)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -633,8 +633,10 @@ def test_equilibrium_form(self, model):
)
assert value(model.rblock[1].equilibrium_constraint["e2"].body) == value(
model.rblock[1].log_k_eq["e2"]
- model.sblock[1].log_mole_frac_phase_comp["p2", "c1"] * -5
+ model.sblock[1].log_mole_frac_phase_comp["p2", "c2"] * 6
- (
model.sblock[1].log_mole_frac_phase_comp["p2", "c1"] * -5
+ model.sblock[1].log_mole_frac_phase_comp["p2", "c2"] * 6
)
)

@pytest.mark.unit
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -112,7 +112,7 @@ def appr_to_true_species(b, p, j):
# Next, check for inherent reactions
if b.has_inherent_reactions:
for r in b.params.inherent_reaction_idx:
# Get stoichiometric coeffiicient for inherent reactions
# Get stoichiometric coefficient for inherent reactions
gamma = b.params.inherent_reaction_stoichiometry[r, p, j]

if gamma != 0:
Expand Down Expand Up @@ -153,7 +153,6 @@ def _apparent_species_scaling(b):
sf = iscale.get_scaling_factor(
b.flow_mol_phase_comp_true[p, j], default=1, warning=True
)

iscale.constraint_scaling_transform(
b.appr_to_true_species[p, j], sf, overwrite=False
)
Expand Down
9 changes: 8 additions & 1 deletion idaes/models/unit_models/heat_exchanger_ntu.py
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@
from idaes.core.util.math import smooth_min, smooth_max
from idaes.core.solvers import get_solver
from idaes.core.util.exceptions import InitializationError
import idaes.core.util.scaling as iscale
import idaes.logger as idaeslog

__author__ = "Paul Akula, Andrew Lee"
Expand Down Expand Up @@ -371,7 +372,13 @@ def rule_entu(blk, t):

self.heat_duty_constraint = Constraint(self.flowsheet().time, rule=rule_entu)

# TODO : Add scaling methods
def calculate_scaling_factors(self):
super().calculate_scaling_factors()

for t in self.flowsheet().time:
sf_heat = iscale.get_scaling_factor(self.hot_side.heat[t], default=1, warning=True)
iscale.constraint_scaling_transform(self.energy_balance_constraint[t], sf_heat, overwrite=False)
iscale.constraint_scaling_transform(self.heat_duty_constraint[t], sf_heat, overwrite=False)

def initialize_build(
self,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -251,7 +251,7 @@ def return_expression(b, p, j, T):


class PressureSatSolvent:
# Method for calculating saturation pressure ofsolvents
# Method for calculating saturation pressure of solvents
@staticmethod
def build_parameters(cobj):
cobj.pressure_sat_comp_coeff_1 = Var(
Expand Down
35 changes: 27 additions & 8 deletions idaes/models_extra/column_models/solvent_column.py
Original file line number Diff line number Diff line change
Expand Up @@ -394,7 +394,7 @@ def column_cross_section_area_eqn(blk):
self.vapor_phase.length_domain,
initialize=0.9,
units=(pyunits.m) ** 2 / (pyunits.m) ** 3,
doc="Specific interfacial area",
doc="Packing interfacial area per unit of column volume",
)

# Liquid and vapor holdups
Expand Down Expand Up @@ -492,6 +492,7 @@ def mechanical_equilibrium(blk, t, x):
)

# Mass transfer constraints
# "mass_transfer" is a bad name for this variable since it's written with a mole basis. "material" is better
self.interphase_mass_transfer = Var(
self.flowsheet().time,
self.liquid_phase.length_domain,
Expand Down Expand Up @@ -560,10 +561,10 @@ def vapor_mass_transfer_eqn(blk, t, x, j):
elif j in equilibrium_comp:
return (
pyunits.convert(
blk.vapor_phase.mass_transfer_term[t, x, "Vap", j],
to_units=lunits("amount") / lunits("time") / lunits("length"),
-blk.interphase_mass_transfer[t, x, j],
to_units=vunits("amount") / vunits("time") / vunits("length"),
)
== -blk.interphase_mass_transfer[t, x, j]
== blk.vapor_phase.mass_transfer_term[t, x, "Vap", j]
)
else:
return blk.vapor_phase.mass_transfer_term[t, x, "Vap", j] == 0.0
Expand All @@ -575,7 +576,7 @@ def vapor_mass_transfer_eqn(blk, t, x, j):
self.vapor_phase.length_domain,
initialize=100,
units=lunits("power") / lunits("temperature") / lunits("length"),
doc="Vapor-liquid heat transfer coefficient",
doc="Vapor-liquid heat transfer coefficient multiplied by heat transfer area per unit column length",
)

# Heat transfer
Expand Down Expand Up @@ -688,6 +689,13 @@ def pressure_at_interface(blk, t, x, j):
def calculate_scaling_factors(self):
super().calculate_scaling_factors()

vunits = (
self.config.vapor_phase.property_package.get_metadata().get_derived_units
)
lunits = (
self.config.liquid_phase.property_package.get_metadata().get_derived_units
)

# ---------------------------------------------------------------------
# Scale variables
for (t, x, j), v in self.pressure_equil.items():
Expand Down Expand Up @@ -751,23 +759,34 @@ def calculate_scaling_factors(self):
)

for (t, x, j), v in self.liquid_mass_transfer_eqn.items():
zf = self.vapor_phase.length_domain.next(x)
try:
sf = iscale.get_scaling_factor(
self.interphase_mass_transfer[t, x, j], default=1, warning=False
self.interphase_mass_transfer[t, zf, j], default=1, warning=False
)
except KeyError:
# This implies a non-volatile component
sf = 1
sf = iscale.get_scaling_factor(
self.liquid_phase.mass_transfer_term[t, x, "Liq", j], default=1, warning=True
)
iscale.constraint_scaling_transform(v, sf)

for (t, x, j), v in self.vapor_mass_transfer_eqn.items():
try:
sf = iscale.get_scaling_factor(
self.interphase_mass_transfer[t, x, j], default=1, warning=False
)
# Account for the fact that this equation is written on a vapor unit basis
sf_units = pyunits.convert(
1 / (lunits("amount") / lunits("time") / lunits("length")),
1 / (vunits("amount") / vunits("time") / vunits("length"))
)
sf *= sf_units
except KeyError:
# This implies a non-volatile component
sf = 1
sf = iscale.get_scaling_factor(
self.vapor_phase.mass_transfer_term[t, x, "Vap", j], default=1, warning=True
)
iscale.constraint_scaling_transform(v, sf)

for (t, x), v in self.heat_transfer_eqn1.items():
Expand Down
19 changes: 13 additions & 6 deletions idaes/models_extra/column_models/solvent_condenser.py
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@
class SolventCondenserData(UnitModelBlockData):
"""
Condenser unit for solvent column models using separate property packages
for liquid and vpor phases.
for liquid and vapor phases.

Unit model to condense the vapor from the top of a solvent column.
"""
Expand Down Expand Up @@ -321,11 +321,11 @@ def build(self):
flow_basis = self.liquid_phase[t_init].get_material_flow_basis()
if flow_basis == MaterialFlowBasis.molar:
fb = "flow_mole"
elif flow_basis == MaterialFlowBasis.molar:
fb = "flow_mass"
# elif flow_basis == MaterialFlowBasis.mass:
Copy link
Member

Choose a reason for hiding this comment

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

Commented code

# fb = "flow_mass"
else:
raise ConfigurationError(
f"{self.name} SolventCondenser only supports mass or molar "
f"{self.name} SolventCondenser only supports molar "
f"basis for MaterialFlowBasis."
)

Expand All @@ -349,7 +349,7 @@ def rule_material_balance(blk, t, j):
elif j in self.liquid_phase.component_list:
# Non-volatile component
# No mass transfer term
# Set liquid flowrate to an arbitary small value
# Set liquid flowrate to an arbitrary small value
return (
blk.liquid_phase[t].get_material_flow_terms("Liq", j)
== blk.zero_flow_param
Expand Down Expand Up @@ -450,7 +450,14 @@ def calculate_scaling_factors(self):
),
)
elif j in self.liquid_phase.component_list:
iscale.constraint_scaling_transform(v, value(1 / self.zero_flow_param))
iscale.constraint_scaling_transform(
v,
iscale.get_scaling_factor(
self.liquid_phase[t].get_material_flow_terms("Liq",j),
default=1,
warning=True,
),
)
else:
pass # no need to scale this constraint

Expand Down