Coverage for backend/ahuora-builder/src/ahuora_builder/custom/watertap/pressure_exchanger.py: 86%
211 statements
« prev ^ index » next coverage.py v7.10.7, created at 2026-07-22 05:22 +0000
1from idaes.core import declare_process_block_class
2import idaes.logger as idaeslog
3from pyomo.common.collections import ComponentSet
4from pyomo.core.base.component import Component
5from pyomo.core.expr.visitor import identify_variables
6from pyomo.environ import Constraint, check_optimal_termination, value
7from pyomo.repn import generate_standard_repn
8from watertap.core.solvers import get_solver
9from watertap.core.util.initialization import interval_initializer
10from watertap.unit_models.pressure_exchanger import (
11 PressureExchangeType,
12 PressureExchangerData as WaterTAPPressureExchangerData,
13)
14from idaes.core.util.exceptions import ConfigurationError, InitializationError
17@declare_process_block_class("PressureExchanger")
18class PressureExchangerData(WaterTAPPressureExchangerData):
19 """
20 Ahuora wrapper around WaterTAP's Pressure Exchanger unit model.
21 """
23 def build(self):
24 super().build()
26 def initialize_build(
27 self,
28 state_args=None,
29 routine=None,
30 outlvl=idaeslog.NOTSET,
31 solver=None,
32 optarg=None,
33 ):
34 """
35 Initialize the pressure exchanger with robust inlet flow guesses.
37 WaterTAP's pressure exchanger initializer requires equal inlet
38 volumetric flows before solving the unit. Platform "guess" properties
39 are temporary initializer aids, so a flow guess should not prevent this
40 unit from finding its own equal-flow starting point.
41 """
42 init_log = idaeslog.getInitLogger(self.name, outlvl, tag="properties")
43 solve_log = idaeslog.getSolveLogger(self.name, outlvl, tag="properties")
44 opt = get_solver(solver, optarg)
46 released_guess_states = _release_platform_flow_guesses(self)
47 true_flow_spec_states = _inlet_states_with_platform_flow_specs(self)
49 flags_low_in = self.feed_side.properties_in.initialize(
50 outlvl=outlvl,
51 optarg=optarg,
52 solver=solver,
53 state_args=state_args,
54 hold_state=True,
55 )
56 flags_high_in = self.brine_side.properties_in.initialize(
57 outlvl=outlvl,
58 optarg=optarg,
59 solver=solver,
60 state_args=state_args,
61 hold_state=True,
62 )
64 init_log.info_high("Initialize inlets complete")
66 _seed_equal_inlet_flow_values(
67 self,
68 released_guess_states=released_guess_states,
69 true_flow_spec_states=true_flow_spec_states,
70 )
72 if value(self.feed_side.properties_in[0].pressure) > value(
73 self.brine_side.properties_in[0].pressure
74 ):
75 raise ConfigurationError(
76 "Initializing pressure exchanger failed because "
77 "the feed side inlet has a higher pressure "
78 "than the brine side inlet"
79 )
81 if ( 81 ↛ 91line 81 didn't jump to line 91 because the condition on line 81 was never true
82 abs(
83 value(self.feed_side.properties_in[0].flow_vol)
84 - value(self.brine_side.properties_in[0].flow_vol)
85 )
86 / value(self.brine_side.properties_in[0].flow_vol)
87 > 1e-4
88 and not self.config.has_mixing
89 and not self.config.has_leakage
90 ):
91 raise ConfigurationError(
92 "Initializing pressure exchanger failed because "
93 "the volumetric flow rates are not equal for both inlets "
94 + str(value(self.brine_side.properties_out[0].flow_vol))
95 + ","
96 + str(value(self.feed_side.properties_in[0].flow_vol))
97 )
98 else:
99 self.eq_equal_flow_vol.deactivate()
101 def propagate_state(sb1, sb2):
102 state_dict_1 = sb1.define_state_vars()
103 state_dict_2 = sb2.define_state_vars()
104 for key in state_dict_1.keys():
105 if state_dict_1[key].is_indexed():
106 for index in state_dict_1[key].keys():
107 state_dict_2[key][index].value = state_dict_1[key][index].value
108 else:
109 state_dict_2[key].value = state_dict_1[key].value
111 propagate_state(
112 self.feed_side.properties_in[0],
113 self.feed_side.properties_out[0],
114 )
115 if self.config.pressure_exchange_calculation is PressureExchangeType.efficiency: 115 ↛ 122line 115 didn't jump to line 122 because the condition on line 115 was always true
116 self.feed_side.properties_out[0].pressure = self.feed_side.properties_in[
117 0
118 ].pressure.value + self.efficiency_pressure_exchanger[0].value * (
119 self.brine_side.properties_in[0].pressure.value
120 - self.feed_side.properties_in[0].pressure.value
121 )
122 elif (
123 self.config.pressure_exchange_calculation
124 is PressureExchangeType.high_pressure_difference
125 ):
126 self.feed_side.properties_out[0].pressure = (
127 self.brine_side.properties_in[0].pressure.value
128 - self.high_pressure_difference[0].value
129 )
131 propagate_state(
132 self.brine_side.properties_in[0],
133 self.brine_side.properties_out[0],
134 )
135 self.brine_side.properties_out[0].pressure.value = self.feed_side.properties_in[
136 0
137 ].pressure.value
138 init_log.info_high("Initialize outlets complete")
140 interval_initializer(self)
142 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
143 results = opt.solve(self, tee=slc.tee)
144 init_log.info("Initialization complete: {}".format(idaeslog.condition(results)))
146 self.feed_side.properties_in.release_state(flags_low_in)
147 self.brine_side.properties_in.release_state(flags_high_in)
149 if not check_optimal_termination(results): 149 ↛ 150line 149 didn't jump to line 150 because the condition on line 149 was never true
150 raise InitializationError(f"Unit model {self.name} failed to initialize")
152 self.eq_equal_flow_vol.activate()
154 def diagnose(self) -> list[tuple[Component, str]]:
155 """
156 Return user-facing hints for common pressure exchanger specification issues.
158 The pressure exchanger expects low-pressure feed on ``feed_side`` and
159 high-pressure retentate/brine on ``brine_side``. ``deltaP`` and ``work``
160 are calculated performance values in the default efficiency mode, so
161 fixing them alongside efficiency commonly over-specifies the unit.
162 """
163 problems: list[tuple[Component, str]] = []
165 feed_in = self.feed_side.properties_in[0]
166 feed_out = self.feed_side.properties_out[0]
167 brine_in = self.brine_side.properties_in[0]
168 brine_out = self.brine_side.properties_out[0]
170 feed_in_pressure = _safe_value(feed_in.pressure)
171 feed_out_pressure = _safe_value(feed_out.pressure)
172 brine_in_pressure = _safe_value(brine_in.pressure)
173 brine_out_pressure = _safe_value(brine_out.pressure)
175 if (
176 feed_in_pressure is not None
177 and brine_in_pressure is not None
178 and feed_in_pressure > brine_in_pressure
179 ):
180 problems.append(
181 (
182 feed_in.pressure,
183 "The pressure exchanger feed side inlet pressure is higher "
184 "than the brine side inlet pressure. The feed side should "
185 "receive low-pressure feed, and the brine side should receive "
186 "high-pressure RO retentate. Check that the low-pressure feed "
187 "is connected to feed_inlet and RO retentate is connected to "
188 "brine_inlet.",
189 )
190 )
191 elif (
192 feed_in_pressure is not None
193 and brine_in_pressure is not None
194 and brine_in_pressure - feed_in_pressure <= 1e-6
195 ):
196 problems.append(
197 (
198 brine_in.pressure,
199 "The brine side inlet pressure is not higher than the feed "
200 "side inlet pressure, so there is no pressure difference to "
201 "recover. The brine side should be the high-pressure stream.",
202 )
203 )
205 if (
206 feed_in_pressure is not None
207 and feed_out_pressure is not None
208 and feed_out_pressure <= feed_in_pressure
209 ):
210 problems.append(
211 (
212 feed_out.pressure,
213 "The feed outlet pressure is not higher than the feed inlet "
214 "pressure. A pressure exchanger should raise pressure on the "
215 "feed side.",
216 )
217 )
219 if (
220 brine_in_pressure is not None
221 and brine_out_pressure is not None
222 and brine_out_pressure >= brine_in_pressure
223 ):
224 problems.append(
225 (
226 brine_out.pressure,
227 "The brine outlet pressure is not lower than the brine inlet "
228 "pressure. A pressure exchanger should recover pressure from "
229 "the high-pressure brine side.",
230 )
231 )
233 feed_flow_vol = _safe_value(feed_in.flow_vol)
234 brine_flow_vol = _safe_value(brine_in.flow_vol)
235 if ( 235 ↛ 260line 235 didn't jump to line 260 because the condition on line 235 was always true
236 feed_flow_vol is not None
237 and brine_flow_vol is not None
238 and min(abs(feed_flow_vol), abs(brine_flow_vol)) > 1e-12
239 ):
240 flow_ratio = abs(feed_flow_vol / brine_flow_vol)
241 if flow_ratio > 100:
242 problems.append(
243 (
244 feed_in.flow_vol,
245 "The feed side volumetric flow is more than 100 times "
246 "the brine side flow. With no leakage or mixing, the "
247 "pressure exchanger expects similar volumetric flows.",
248 )
249 )
250 elif flow_ratio < 0.01:
251 problems.append(
252 (
253 brine_in.flow_vol,
254 "The brine side volumetric flow is more than 100 times "
255 "the feed side flow. With no leakage or mixing, the "
256 "pressure exchanger expects similar volumetric flows.",
257 )
258 )
260 efficiency = _safe_value(self.efficiency_pressure_exchanger[0])
261 if efficiency is not None: 261 ↛ 280line 261 didn't jump to line 280 because the condition on line 261 was always true
262 if efficiency <= 0 or efficiency > 1:
263 problems.append(
264 (
265 self.efficiency_pressure_exchanger[0],
266 "Pressure exchanger efficiency must be greater than 0 "
267 "and no more than 1.",
268 )
269 )
270 elif efficiency < 0.5:
271 problems.append(
272 (
273 self.efficiency_pressure_exchanger[0],
274 "Pressure exchanger efficiency is very low. Confirm this "
275 "is intentional; typical pressure exchangers are much "
276 "closer to 1.",
277 )
278 )
280 if not _is_specified(self, self.efficiency_pressure_exchanger[0]) and not hasattr(
281 self, "high_pressure_difference"
282 ):
283 problems.append(
284 (
285 self.efficiency_pressure_exchanger[0],
286 "Pressure exchanger efficiency is not specified. For a simulation "
287 "case, specify efficiency or use the high-pressure-difference "
288 "calculation mode with its pressure differences specified.",
289 )
290 )
292 for component, label in (
293 (self.feed_side.deltaP[0], "feed side pressure change"),
294 (self.brine_side.deltaP[0], "brine side pressure change"),
295 ):
296 if component.fixed:
297 problems.append(
298 (
299 component,
300 f"The {label} is fixed, but it is calculated from the "
301 "pressure exchanger equations. Fixing it alongside "
302 "efficiency or outlet pressure can over-specify the unit.",
303 )
304 )
306 if feed_out.pressure.fixed and brine_out.pressure.fixed:
307 problems.append(
308 (
309 feed_out.pressure,
310 "Both pressure exchanger outlet pressures are fixed. Outlet "
311 "pressures are coupled to the pressure transfer equations, so "
312 "fixing both can over-specify the unit.",
313 )
314 )
316 for constraint in _constraints_fixing_expression(self, self.feed_side.work[0]):
317 problems.append(
318 (
319 constraint,
320 "Feed side mechanical work is a calculated expression. A "
321 "constraint fixing it can over-specify the pressure exchanger.",
322 )
323 )
325 for constraint in _constraints_fixing_expression(self, self.brine_side.work[0]):
326 problems.append(
327 (
328 constraint,
329 "Brine side mechanical work is a calculated expression. A "
330 "constraint fixing it can over-specify the pressure exchanger.",
331 )
332 )
334 return problems
337def _safe_value(component) -> float | None:
338 try:
339 return value(component, exception=False)
340 except Exception:
341 return None
344def _constraints_fixing_expression(block, expression):
345 expr_id = id(expression)
346 for constraint in block.parent_block().component_data_objects(
347 Constraint, active=True, descend_into=True
348 ):
349 if id(constraint.body) == expr_id:
350 yield constraint
353def _is_specified(block, component) -> bool:
354 return component.fixed or _has_single_variable_equality_constraint(block, component)
357def _has_single_variable_equality_constraint(block, component) -> bool:
358 component_id = id(component)
359 for constraint in block.parent_block().component_data_objects(
360 Constraint, active=True, descend_into=True
361 ):
362 if _is_descendant_of(constraint.parent_block(), block):
363 continue
364 if not constraint.equality: 364 ↛ 365line 364 didn't jump to line 365 because the condition on line 364 was never true
365 continue
367 repn = generate_standard_repn(constraint.body)
368 if ( 368 ↛ 375line 368 didn't jump to line 375 because the condition on line 368 was always true
369 repn.is_linear()
370 and len(repn.linear_vars) == 1
371 and id(repn.linear_vars[0]) == component_id
372 ):
373 return True
375 variables = list(identify_variables(constraint.body, include_fixed=False))
376 if len(variables) == 1 and id(variables[0]) == component_id:
377 return True
379 return False
382def _is_descendant_of(child_block, parent_block) -> bool:
383 block = child_block
384 while block is not None:
385 if block is parent_block:
386 return True
387 block = block.parent_block()
388 return False
391def _release_platform_flow_guesses(unit) -> ComponentSet:
392 """
393 Unfix platform flow guesses on inlet state blocks before WaterTAP initializes.
395 Flexible state blocks normally reactivate values registered through
396 ``constrain_component`` during their own initialization. That is right for
397 real specifications, but platform guesses should remain free after their
398 value has seeded the model.
399 """
400 guess_components = _platform_guess_components(unit)
401 released_states = ComponentSet()
402 if not guess_components:
403 return released_states
405 for state_block in _pressure_exchanger_inlet_states(unit):
406 vars_to_deactivate = state_block.__dict__.get("vars_to_deactivate")
407 if vars_to_deactivate is None: 407 ↛ 408line 407 didn't jump to line 408 because the condition on line 407 was never true
408 continue
410 retained_vars = []
411 for component in vars_to_deactivate:
412 if component in guess_components and _is_flow_component(component):
413 component.unfix()
414 released_states.add(state_block)
415 else:
416 retained_vars.append(component)
417 state_block.__dict__["vars_to_deactivate"] = retained_vars
419 for component in guess_components:
420 if _is_flow_component(component) and hasattr(component, "unfix"): 420 ↛ 419line 420 didn't jump to line 419 because the condition on line 420 was always true
421 component.unfix()
423 return released_states
426def _seed_equal_inlet_flow_values(
427 unit,
428 *,
429 released_guess_states: ComponentSet,
430 true_flow_spec_states: ComponentSet,
431) -> None:
432 if unit.config.has_mixing or unit.config.has_leakage: 432 ↛ 433line 432 didn't jump to line 433 because the condition on line 432 was never true
433 return
435 feed = unit.feed_side.properties_in[0]
436 brine = unit.brine_side.properties_in[0]
437 feed_flow_vol = _safe_value(feed.flow_vol)
438 brine_flow_vol = _safe_value(brine.flow_vol)
439 if ( 439 ↛ 445line 439 didn't jump to line 445 because the condition on line 439 was never true
440 feed_flow_vol is None
441 or brine_flow_vol is None
442 or abs(feed_flow_vol) <= 1e-12
443 or abs(brine_flow_vol) <= 1e-12
444 ):
445 return
447 if abs(feed_flow_vol - brine_flow_vol) / abs(brine_flow_vol) <= 1e-4:
448 return
450 feed_can_move = feed in released_guess_states or feed not in true_flow_spec_states
451 brine_can_move = brine in released_guess_states or brine not in true_flow_spec_states
452 if feed_can_move: 452 ↛ 454line 452 didn't jump to line 454 because the condition on line 452 was always true
453 _scale_flow_state_to_volume(feed, brine_flow_vol)
454 elif brine_can_move:
455 _scale_flow_state_to_volume(brine, feed_flow_vol)
458def _pressure_exchanger_inlet_states(unit):
459 return (
460 unit.feed_side.properties_in[0],
461 unit.brine_side.properties_in[0],
462 )
465def _platform_guess_components(unit) -> ComponentSet:
466 try:
467 guess_groups = unit.flowsheet().guess_vars
468 except AttributeError:
469 return ComponentSet()
471 components = ComponentSet()
472 inlet_states = ComponentSet(_pressure_exchanger_inlet_states(unit))
473 for guess_group in guess_groups:
474 for component in guess_group:
475 if component.parent_block() in inlet_states: 475 ↛ 474line 475 didn't jump to line 474 because the condition on line 475 was always true
476 components.add(component)
477 return components
480def _inlet_states_with_platform_flow_specs(unit) -> ComponentSet:
481 states = ComponentSet()
482 for state_block in _pressure_exchanger_inlet_states(unit):
483 for component in state_block.__dict__.get("vars_to_deactivate", []):
484 if _is_flow_component(component):
485 states.add(state_block)
486 break
487 return states
490def _is_flow_component(component) -> bool:
491 return "flow" in component.local_name
494def _scale_flow_state_to_volume(state_block, target_flow_vol: float) -> None:
495 current_flow_vol = _safe_value(state_block.flow_vol)
496 if current_flow_vol is None or abs(current_flow_vol) <= 1e-12: 496 ↛ 497line 496 didn't jump to line 497 because the condition on line 496 was never true
497 return
499 flow_mass_vars = list(state_block.flow_mass_phase_comp.values())
500 current_total_flow_mass = sum(
501 _safe_value(component) or 0 for component in flow_mass_vars
502 )
503 target_total_flow_mass = current_total_flow_mass * target_flow_vol / current_flow_vol
504 if current_total_flow_mass <= 1e-12: 504 ↛ 505line 504 didn't jump to line 505 because the condition on line 504 was never true
505 per_component_flow_mass = target_total_flow_mass / len(flow_mass_vars)
506 for component in flow_mass_vars:
507 component.set_value(per_component_flow_mass)
508 return
510 scale_factor = target_total_flow_mass / current_total_flow_mass
511 for component in flow_mass_vars:
512 component.set_value((_safe_value(component) or 0) * scale_factor)
514 if hasattr(state_block, "flow_vol_phase"): 514 ↛ exitline 514 didn't return from function '_scale_flow_state_to_volume' because the condition on line 514 was always true
515 for component in state_block.flow_vol_phase.values():
516 component.set_value(target_flow_vol / len(state_block.flow_vol_phase))