diff --git a/file24.py b/file24.py new file mode 100644 index 0000000..af80553 --- /dev/null +++ b/file24.py @@ -0,0 +1,98 @@ +fs.compressor.control_volume.properties_in[0.0].flow_mol: 1 fixed +fs.compressor.control_volume.properties_in[0.0].mole_frac_comp[oxygen]: 0.5 fixed +fs.compressor.control_volume.properties_in[0.0].mole_frac_comp[nitrogen]: 0.5 fixed +fs.compressor.control_volume.properties_in[0.0].temperature: 298.08186264499585 unfixed +fs.compressor.control_volume.properties_in[0.0].pressure: 100000 fixed +fs.compressor.control_volume.properties_out[0.0].flow_mol: 6.57562972162535e-10 unfixed +fs.compressor.control_volume.properties_out[0.0].mole_frac_comp[oxygen]: 0.6267749501040214 unfixed +fs.compressor.control_volume.properties_out[0.0].mole_frac_comp[nitrogen]: 0.37322504989597827 unfixed +fs.compressor.control_volume.properties_out[0.0].temperature: 540.815233487118 unfixed +fs.compressor.control_volume.properties_out[0.0].pressure: 52495.189577299374 unfixed +fs.compressor.control_volume.work[0.0]: 12.500000194876284 unfixed +fs.compressor.control_volume.deltaP[0.0]: -47504.810422700626 unfixed +fs.compressor.ratioP[0.0]: 0.5249518957729937 unfixed +fs.compressor.efficiency_isentropic[0.0]: 0.8 fixed +fs.compressor.work_isentropic[0.0]: 10.000000155903408 unfixed +fs.compressor.control_volume.properties_in[0.0].flow_mol_phase[Liq]: 0.023005058626483423 unfixed +fs.compressor.control_volume.properties_in[0.0].flow_mol_phase[Vap]: 0.9769949413735161 unfixed +fs.compressor.control_volume.properties_in[0.0].mole_frac_phase_comp[Liq,oxygen]: 0.7814553536292846 unfixed +fs.compressor.control_volume.properties_in[0.0].mole_frac_phase_comp[Liq,nitrogen]: 0.18913795627561378 unfixed +fs.compressor.control_volume.properties_in[0.0].mole_frac_phase_comp[Vap,oxygen]: 0.4854838496927152 unfixed +fs.compressor.control_volume.properties_in[0.0].mole_frac_phase_comp[Vap,nitrogen]: 0.4851094602121835 unfixed +fs.compressor.control_volume.properties_in[0.0].phase_frac[Liq]: 0.023005058631352215 unfixed +fs.compressor.control_volume.properties_in[0.0].phase_frac[Vap]: 0.9769949413781533 unfixed +fs.compressor.control_volume.properties_in[0.0]._teq[Vap,Liq]: 85.35128696195348 unfixed +fs.compressor.control_volume.properties_in[0.0]._t1_Vap_Liq: 298.0818627602398 unfixed +fs.compressor.control_volume.properties_in[0.0].temperature_bubble[Vap,Liq]: 81.15070769322699 unfixed +fs.compressor.control_volume.properties_in[0.0]._mole_frac_tbub[Vap,Liq,oxygen]: 0.1963099228994071 unfixed +fs.compressor.control_volume.properties_in[0.0]._mole_frac_tbub[Vap,Liq,nitrogen]: 0.8036900771005928 unfixed +fs.compressor.control_volume.properties_in[0.0].log_mole_frac_comp[oxygen]: -0.6931471805599475 unfixed +fs.compressor.control_volume.properties_in[0.0].log_mole_frac_comp[nitrogen]: -0.6931471805599458 unfixed +fs.compressor.control_volume.properties_in[0.0].log_mole_frac_tbub[Vap,Liq,oxygen]: -1.6280606293114597 unfixed +fs.compressor.control_volume.properties_in[0.0].log_mole_frac_tbub[Vap,Liq,nitrogen]: -0.21854156036127984 unfixed +fs.compressor.control_volume.properties_in[0.0].temperature_dew[Vap,Liq]: 85.35128696224727 unfixed +fs.compressor.control_volume.properties_in[0.0]._mole_frac_tdew[Vap,Liq,oxygen]: 0.8080791333675067 unfixed +fs.compressor.control_volume.properties_in[0.0]._mole_frac_tdew[Vap,Liq,nitrogen]: 0.1919208666324939 unfixed +fs.compressor.control_volume.properties_in[0.0].log_mole_frac_tdew[Vap,Liq,oxygen]: -0.2130952879205655 unfixed +fs.compressor.control_volume.properties_in[0.0].log_mole_frac_tdew[Vap,Liq,nitrogen]: -1.6506721448684876 unfixed +fs.compressor.control_volume.properties_in[0.0].log_mole_frac_phase_comp[Liq,oxygen]: -0.24659725982574596 unfixed +fs.compressor.control_volume.properties_in[0.0].log_mole_frac_phase_comp[Liq,nitrogen]: -1.6652786028057271 unfixed +fs.compressor.control_volume.properties_in[0.0].log_mole_frac_phase_comp[Vap,oxygen]: -0.7226092571125734 unfixed +fs.compressor.control_volume.properties_in[0.0].log_mole_frac_phase_comp[Vap,nitrogen]: -0.7233807223499293 unfixed +fs.compressor.control_volume.properties_out[0.0].flow_mol_phase[Liq]: 0.0 unfixed +fs.compressor.control_volume.properties_out[0.0].flow_mol_phase[Vap]: 9.559685150849945e-09 unfixed +fs.compressor.control_volume.properties_out[0.0].mole_frac_phase_comp[Liq,oxygen]: 3.436473043731984e-07 unfixed +fs.compressor.control_volume.properties_out[0.0].mole_frac_phase_comp[Liq,nitrogen]: 0.021785521327024382 unfixed +fs.compressor.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,oxygen]: 0.02149969775415102 unfixed +fs.compressor.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,nitrogen]: 0.00028616728882837814 unfixed +fs.compressor.control_volume.properties_out[0.0].phase_frac[Liq]: 4.784612919644563e-06 unfixed +fs.compressor.control_volume.properties_out[0.0].phase_frac[Vap]: 14.719623584013826 unfixed +fs.compressor.control_volume.properties_out[0.0]._teq[Vap,Liq]: 176.86836842529004 unfixed +fs.compressor.control_volume.properties_out[0.0]._t1_Vap_Liq: 540.8152335410161 unfixed +fs.compressor.control_volume.properties_out[0.0].temperature_bubble[Vap,Liq]: 76.9776019724163 unfixed +fs.compressor.control_volume.properties_out[0.0]._mole_frac_tbub[Vap,Liq,oxygen]: 0.2548973116378595 unfixed +fs.compressor.control_volume.properties_out[0.0]._mole_frac_tbub[Vap,Liq,nitrogen]: 0.7451026883621408 unfixed +fs.compressor.control_volume.properties_out[0.0].log_mole_frac_comp[oxygen]: -0.46716773403401146 unfixed +fs.compressor.control_volume.properties_out[0.0].log_mole_frac_comp[nitrogen]: -0.9855736903537607 unfixed +fs.compressor.control_volume.properties_out[0.0].log_mole_frac_tbub[Vap,Liq,oxygen]: -1.3668945143882842 unfixed +fs.compressor.control_volume.properties_out[0.0].log_mole_frac_tbub[Vap,Liq,nitrogen]: -0.29423323337348073 unfixed +fs.compressor.control_volume.properties_out[0.0].temperature_dew[Vap,Liq]: 176.86836842546165 unfixed +fs.compressor.control_volume.properties_out[0.0]._mole_frac_tdew[Vap,Liq,oxygen]: 0.6267749501040248 unfixed +fs.compressor.control_volume.properties_out[0.0]._mole_frac_tdew[Vap,Liq,nitrogen]: 0.3732250498959747 unfixed +fs.compressor.control_volume.properties_out[0.0].log_mole_frac_tdew[Vap,Liq,oxygen]: -0.4671677340340085 unfixed +fs.compressor.control_volume.properties_out[0.0].log_mole_frac_tdew[Vap,Liq,nitrogen]: -0.985573690353764 unfixed +fs.compressor.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Liq,oxygen]: -14.883739171323969 unfixed +fs.compressor.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Liq,nitrogen]: -3.826509689029257 unfixed +fs.compressor.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Vap,oxygen]: -3.8397164016317236 unfixed +fs.compressor.control_volume.properties_out[0.0].log_mole_frac_phase_comp[Vap,nitrogen]: -8.158933965396713 unfixed +fs.compressor.properties_isentropic[0.0].flow_mol: 3.4799210529450807e-09 unfixed +fs.compressor.properties_isentropic[0.0].mole_frac_comp[oxygen]: 0.506165142531512 unfixed +fs.compressor.properties_isentropic[0.0].mole_frac_comp[nitrogen]: 0.4938348574684874 unfixed +fs.compressor.properties_isentropic[0.0].pressure: 52495.189577299374 unfixed +fs.compressor.properties_isentropic[0.0].temperature: 300 unfixed +fs.compressor.properties_isentropic[0.0].flow_mol_phase[Liq]: 8.565243957053452e-11 unfixed +fs.compressor.properties_isentropic[0.0].flow_mol_phase[Vap]: 3.1303204846938315e-09 unfixed +fs.compressor.properties_isentropic[0.0].mole_frac_phase_comp[Liq,oxygen]: 0.8312971777273821 unfixed +fs.compressor.properties_isentropic[0.0].mole_frac_phase_comp[Liq,nitrogen]: 0.16404817682570494 unfixed +fs.compressor.properties_isentropic[0.0].mole_frac_phase_comp[Vap,oxygen]: 0.5038115958529246 unfixed +fs.compressor.properties_isentropic[0.0].mole_frac_phase_comp[Vap,nitrogen]: 0.4915337587001636 unfixed +fs.compressor.properties_isentropic[0.0].phase_frac[Liq]: 0.044004181848622224 unfixed +fs.compressor.properties_isentropic[0.0].phase_frac[Vap]: 1.0013268160053315 unfixed +fs.compressor.properties_isentropic[0.0]._teq[Vap,Liq]: 80.08261235742805 unfixed +fs.compressor.properties_isentropic[0.0]._t1_Vap_Liq: 300.0000001114594 unfixed +fs.compressor.properties_isentropic[0.0].temperature_bubble[Vap,Liq]: 75.7032281335816 unfixed +fs.compressor.properties_isentropic[0.0]._mole_frac_tbub[Vap,Liq,oxygen]: 0.17517382967610917 unfixed +fs.compressor.properties_isentropic[0.0]._mole_frac_tbub[Vap,Liq,nitrogen]: 0.8248261703238914 unfixed +fs.compressor.properties_isentropic[0.0].log_mole_frac_comp[oxygen]: -0.6808922943019795 unfixed +fs.compressor.properties_isentropic[0.0].log_mole_frac_comp[nitrogen]: -0.7055541143062826 unfixed +fs.compressor.properties_isentropic[0.0].log_mole_frac_tbub[Vap,Liq,oxygen]: -1.7419764856320588 unfixed +fs.compressor.properties_isentropic[0.0].log_mole_frac_tbub[Vap,Liq,nitrogen]: -0.192582617486082 unfixed +fs.compressor.properties_isentropic[0.0].temperature_dew[Vap,Liq]: 80.08261235771225 unfixed +fs.compressor.properties_isentropic[0.0]._mole_frac_tdew[Vap,Liq,oxygen]: 0.8355815665536349 unfixed +fs.compressor.properties_isentropic[0.0]._mole_frac_tdew[Vap,Liq,nitrogen]: 0.16441843344636473 unfixed +fs.compressor.properties_isentropic[0.0].log_mole_frac_tdew[Vap,Liq,oxygen]: -0.17962730967513357 unfixed +fs.compressor.properties_isentropic[0.0].log_mole_frac_tdew[Vap,Liq,nitrogen]: -1.8053406770595029 unfixed +fs.compressor.properties_isentropic[0.0].log_mole_frac_phase_comp[Liq,oxygen]: -0.1847679334618234 unfixed +fs.compressor.properties_isentropic[0.0].log_mole_frac_phase_comp[Liq,nitrogen]: -1.8075951331649516 unfixed +fs.compressor.properties_isentropic[0.0].log_mole_frac_phase_comp[Vap,oxygen]: -0.6855528985501691 unfixed +fs.compressor.properties_isentropic[0.0].log_mole_frac_phase_comp[Vap,nitrogen]: -0.7102246567075116 unfixed diff --git a/file_fs.state.txt b/file_fs.state.txt index 12cafe5..8a8723f 100644 --- a/file_fs.state.txt +++ b/file_fs.state.txt @@ -1,42 +1,32 @@ fs.state[0].flow_mol: 1 fixed -fs.state[0].mole_frac_comp[oxygen]: 0.33 fixed -fs.state[0].mole_frac_comp[argon]: 0.33 fixed -fs.state[0].mole_frac_comp[nitrogen]: 0.33 fixed -fs.state[0].pressure: 200000 fixed -fs.state[0].enth_mol: 4000 fixed -fs.state[0].flow_mol_phase[Liq]: 0.05503450613138045 unfixed -fs.state[0].flow_mol_phase[Vap]: 0.9449654938686196 unfixed -fs.state[0].mole_frac_phase_comp[Liq,oxygen]: 0.47077791836720434 unfixed -fs.state[0].mole_frac_phase_comp[Liq,argon]: 0.36891049697959266 unfixed -fs.state[0].mole_frac_phase_comp[Liq,nitrogen]: 0.15031158465320293 unfixed -fs.state[0].mole_frac_phase_comp[Vap,oxygen]: 0.3218010422538228 unfixed -fs.state[0].mole_frac_phase_comp[Vap,argon]: 0.32773388379929397 unfixed -fs.state[0].mole_frac_phase_comp[Vap,nitrogen]: 0.3404650739468832 unfixed -fs.state[0].temperature: 92.4653653053216 unfixed -fs.state[0].phase_frac[Liq]: 0.05503450613138045 unfixed -fs.state[0].phase_frac[Vap]: 0.9449654938686196 unfixed -fs.state[0]._teq[Vap,Liq]: 92.46512005288838 unfixed -fs.state[0]._t1_Vap_Liq: 92.46537475236082 unfixed -fs.state[0].temperature_bubble[Vap,Liq]: 89.81904298694354 unfixed -fs.state[0]._mole_frac_tbub[Vap,Liq,oxygen]: 0.17962142058398192 unfixed -fs.state[0]._mole_frac_tbub[Vap,Liq,argon]: 0.23129754896474208 unfixed -fs.state[0]._mole_frac_tbub[Vap,Liq,nitrogen]: 0.589081030451276 unfixed -fs.state[0].log_mole_frac_comp[oxygen]: -1.1086626245216111 unfixed -fs.state[0].log_mole_frac_comp[argon]: -1.1086626245216111 unfixed -fs.state[0].log_mole_frac_comp[nitrogen]: -1.1086626245216111 unfixed -fs.state[0].log_mole_frac_tbub[Vap,Liq,oxygen]: -1.716903861940957 unfixed -fs.state[0].log_mole_frac_tbub[Vap,Liq,argon]: -1.4640503065810795 unfixed -fs.state[0].log_mole_frac_tbub[Vap,Liq,nitrogen]: -0.529191531870467 unfixed -fs.state[0].temperature_dew[Vap,Liq]: 92.46536544012613 unfixed -fs.state[0]._mole_frac_tdew[Vap,Liq,oxygen]: 0.4839403861534373 unfixed -fs.state[0]._mole_frac_tdew[Vap,Liq,argon]: 0.371437904909194 unfixed -fs.state[0]._mole_frac_tdew[Vap,Liq,nitrogen]: 0.14462170893736873 unfixed -fs.state[0].log_mole_frac_tdew[Vap,Liq,oxygen]: -0.7257935489559108 unfixed -fs.state[0].log_mole_frac_tdew[Vap,Liq,argon]: -0.9903735757419188 unfixed -fs.state[0].log_mole_frac_tdew[Vap,Liq,nitrogen]: -1.9336338495621757 unfixed -fs.state[0].log_mole_frac_phase_comp[Liq,oxygen]: -0.7533688070595927 unfixed -fs.state[0].log_mole_frac_phase_comp[Liq,argon]: -0.9972012199738397 unfixed -fs.state[0].log_mole_frac_phase_comp[Liq,nitrogen]: -1.8950449083258598 unfixed -fs.state[0].log_mole_frac_phase_comp[Vap,oxygen]: -1.1338218056042926 unfixed -fs.state[0].log_mole_frac_phase_comp[Vap,argon]: -1.1155533297844114 unfixed -fs.state[0].log_mole_frac_phase_comp[Vap,nitrogen]: -1.077442731496745 unfixed +fs.state[0].mole_frac_comp[oxygen]: 0.5 fixed +fs.state[0].mole_frac_comp[nitrogen]: 0.5 fixed +fs.state[0].pressure: 100000 fixed +fs.state[0].enth_mol: -10 fixed +fs.state[0].flow_mol_phase[Liq]: 7.071166685008625e-05 unfixed +fs.state[0].flow_mol_phase[Vap]: 0.9999292883331499 unfixed +fs.state[0].mole_frac_phase_comp[Liq,oxygen]: 0.8080656169942638 unfixed +fs.state[0].mole_frac_phase_comp[Liq,nitrogen]: 0.19193438300573623 unfixed +fs.state[0].mole_frac_phase_comp[Vap,oxygen]: 0.49997821518309227 unfixed +fs.state[0].mole_frac_phase_comp[Vap,nitrogen]: 0.5000217848169077 unfixed +fs.state[0].temperature: 85.3512853489739 unfixed +fs.state[0].phase_frac[Liq]: 7.071166685008625e-05 unfixed +fs.state[0].phase_frac[Vap]: 0.9999292883331499 unfixed +fs.state[0]._teq[Vap,Liq]: 85.35103912197744 unfixed +fs.state[0]._t1_Vap_Liq: 85.35129130052786 unfixed +fs.state[0].temperature_bubble[Vap,Liq]: 81.15070769322693 unfixed +fs.state[0]._mole_frac_tbub[Vap,Liq,oxygen]: 0.19630992289940427 unfixed +fs.state[0]._mole_frac_tbub[Vap,Liq,nitrogen]: 0.8036900771005957 unfixed +fs.state[0].log_mole_frac_comp[oxygen]: -0.6931471805599453 unfixed +fs.state[0].log_mole_frac_comp[nitrogen]: -0.6931471805599453 unfixed +fs.state[0].log_mole_frac_tbub[Vap,Liq,oxygen]: -1.6280606293114737 unfixed +fs.state[0].log_mole_frac_tbub[Vap,Liq,nitrogen]: -0.21854156036127612 unfixed +fs.state[0].temperature_dew[Vap,Liq]: 85.35128696224731 unfixed +fs.state[0]._mole_frac_tdew[Vap,Liq,oxygen]: 0.808079133367472 unfixed +fs.state[0]._mole_frac_tdew[Vap,Liq,nitrogen]: 0.19192086663252803 unfixed +fs.state[0].log_mole_frac_tdew[Vap,Liq,oxygen]: -0.21309528792060936 unfixed +fs.state[0].log_mole_frac_tdew[Vap,Liq,nitrogen]: -1.6506721448683128 unfixed +fs.state[0].log_mole_frac_phase_comp[Liq,oxygen]: -0.21311201460700968 unfixed +fs.state[0].log_mole_frac_phase_comp[Liq,nitrogen]: -1.6506017205444223 unfixed +fs.state[0].log_mole_frac_phase_comp[Vap,oxygen]: -0.6931907511429448 unfixed +fs.state[0].log_mole_frac_phase_comp[Vap,nitrogen]: -0.6931036118752588 unfixed diff --git a/property_packages/modular/modular_custom_utils.py b/property_packages/modular/modular_custom_utils.py deleted file mode 100644 index 7452bae..0000000 --- a/property_packages/modular/modular_custom_utils.py +++ /dev/null @@ -1,399 +0,0 @@ -from pyomo.environ import ( - Block, - check_optimal_termination, - Constraint, - exp, - Expression, - log, - Set, - Param, - value, - Var, - units as pyunits, - Reference, - SolverFactory -) -from pyomo.common.config import ConfigBlock, ConfigDict, ConfigValue, In, Bool -from pyomo.util.calc_var_value import calculate_variable_from_constraint - -# Import IDAES cores -from idaes.core import ( - declare_process_block_class, - PhysicalParameterBlock, - StateBlockData, - StateBlock, - MaterialFlowBasis, - ElectrolytePropertySet, -) -from idaes.core.base.components import Component, __all_components__ -from idaes.core.base.phases import ( - Phase, - AqueousPhase, - LiquidPhase, - VaporPhase, - __all_phases__, -) -from idaes.core.util.initialization import ( - fix_state_vars, - revert_state_vars, - solve_indexed_blocks, -) -from idaes.core.util.model_statistics import ( - degrees_of_freedom, - number_activated_constraints, -) -from idaes.core.util.exceptions import ( - BurntToast, - ConfigurationError, - PropertyPackageError, - PropertyNotSupportedError, - InitializationError, -) -from idaes.core.util.misc import add_object_reference -from idaes.core.solvers import get_solver -import idaes.logger as idaeslog -import idaes.core.util.scaling as iscale -from idaes.core.initialization.initializer_base import InitializerBase - -from idaes.models.properties.modular_properties.base.generic_reaction import ( - equil_rxn_config, -) -from idaes.models.properties.modular_properties.base.utility import ( - get_method, - get_phase_method, - GenericPropertyPackageError, - StateIndex, - identify_VL_component_list, - # estimate_Tbub, - # estimate_Tdew, - estimate_Pbub, - estimate_Pdew, -) -from idaes.models.properties.modular_properties.phase_equil.bubble_dew import ( - LogBubbleDew, -) -from idaes.models.properties.modular_properties.phase_equil.henry import HenryType - - -def _init_Tbub(blk, T_units): - for pp in blk.params._pe_pairs: - l_phase, _, raoult_comps, henry_comps, _, _ = identify_VL_component_list( - blk, pp - ) - - if raoult_comps == []: - continue - - Tbub0 = _custom_estimate_Tbub(blk, T_units, raoult_comps, henry_comps, l_phase) - - blk.temperature_bubble[pp].set_value(Tbub0) - - for j in raoult_comps: - blk._mole_frac_tbub[pp, j].value = value( - blk.mole_frac_comp[j] - * get_method(blk, "pressure_sat_comp", j)( - blk, blk.params.get_component(j), Tbub0 * T_units - ) - / blk.pressure - ) - if blk.is_property_constructed("log_mole_frac_tbub"): - blk.log_mole_frac_tbub[pp, j].value = value( - log(blk._mole_frac_tbub[pp, j]) - ) - - for j in henry_comps: - blk._mole_frac_tbub[pp, j].value = value( - blk.mole_frac_comp[j] - * blk.params.get_component(j) - .config.henry_component[l_phase]["method"] - .return_expression(blk, l_phase, j, Tbub0 * T_units) - / blk.pressure - ) - if blk.is_property_constructed("log_mole_frac_tbub"): - blk.log_mole_frac_tbub[pp, j].value = value( - log(blk._mole_frac_tbub[pp, j]) - ) - - -def _init_Tdew(blk, T_units): - for pp in blk.params._pe_pairs: - l_phase, _, raoult_comps, henry_comps, _, _ = identify_VL_component_list( - blk, pp - ) - - if raoult_comps == []: - continue - - Tdew0 = _custom_estimate_Tdew(blk, T_units, raoult_comps, henry_comps, l_phase) - - blk.temperature_dew[pp].set_value(Tdew0) - - for j in raoult_comps: - blk._mole_frac_tdew[pp, j].value = value( - blk.mole_frac_comp[j] - * blk.pressure - / get_method(blk, "pressure_sat_comp", j)( - blk, blk.params.get_component(j), Tdew0 * T_units - ) - ) - if blk.is_property_constructed("log_mole_frac_tdew"): - blk.log_mole_frac_tdew[pp, j].value = value( - log(blk._mole_frac_tdew[pp, j]) - ) - for j in henry_comps: - blk._mole_frac_tdew[pp, j].value = value( - blk.mole_frac_comp[j] - * blk.pressure - / blk.params.get_component(j) - .config.henry_component[l_phase]["method"] - .return_expression(blk, l_phase, j, Tdew0 * T_units) - ) - if blk.is_property_constructed("log_mole_frac_tdew"): - blk.log_mole_frac_tdew[pp, j].value = value( - log(blk._mole_frac_tdew[pp, j]) - ) - - -def _init_Pbub(blk): - for pp in blk.params._pe_pairs: - l_phase, _, raoult_comps, henry_comps, _, _ = identify_VL_component_list( - blk, pp - ) - - if raoult_comps == []: - continue - - blk.pressure_bubble[pp].set_value( - estimate_Pbub(blk, raoult_comps, henry_comps, l_phase) - ) - - for j in raoult_comps: - blk._mole_frac_pbub[pp, j].value = value( - blk.mole_frac_comp[j] - * blk.pressure_sat_comp[j] - / blk.pressure_bubble[pp] - ) - if blk.is_property_constructed("log_mole_frac_pbub"): - blk.log_mole_frac_pbub[pp, j].value = value( - log(blk._mole_frac_pbub[pp, j]) - ) - for j in henry_comps: - blk._mole_frac_pbub[pp, j].value = value( - blk.mole_frac_comp[j] * blk.henry[l_phase, j] / blk.pressure_bubble[pp] - ) - if blk.is_property_constructed("log_mole_frac_pbub"): - blk.log_mole_frac_pbub[pp, j].value = value( - log(blk._mole_frac_pbub[pp, j]) - ) - - -def _init_Pdew(blk): - for pp in blk.params._pe_pairs: - l_phase, _, raoult_comps, henry_comps, _, _ = identify_VL_component_list( - blk, pp - ) - - if raoult_comps == []: - continue - - blk.pressure_dew[pp].set_value( - estimate_Pdew(blk, raoult_comps, henry_comps, l_phase) - ) - - for j in raoult_comps: - blk._mole_frac_pdew[pp, j].value = value( - blk.mole_frac_comp[j] * blk.pressure_dew[pp] / blk.pressure_sat_comp[j] - ) - if blk.is_property_constructed("log_mole_frac_pdew"): - blk.log_mole_frac_pdew[pp, j].value = value( - log(blk._mole_frac_pdew[pp, j]) - ) - for j in henry_comps: - blk._mole_frac_pdew[pp, j].value = value( - blk.mole_frac_comp[j] * blk.pressure_dew[pp] / blk.henry[l_phase, j] - ) - if blk.is_property_constructed("log_mole_frac_pdew"): - blk.log_mole_frac_pdew[pp, j].value = value( - log(blk._mole_frac_pdew[pp, j]) - ) - - - -from enum import Enum - -from pyomo.environ import units as pyunits, value - -from idaes.core.util.exceptions import ( - BurntToast, - ConfigurationError, - PropertyPackageError, -) -import idaes.logger as idaeslog - -# Set up logger -_log = idaeslog.getLogger(__name__) - - -TOL = 1e-1 -MAX_ITER = 30 - - -def _custom_estimate_Tbub(blk, T_units, raoult_comps, henry_comps, liquid_phase): - """ - Function to estimate bubble point temperature - - Args: - blk: StateBlock to use - T_units: units of temperature - raoult_comps: list of components that follow Raoult's Law - henry_comps: list of components that follow Henry's Law - liquid_phase: name of liquid phase - - Returns: - Estimated bubble point temperature as a float. - - """ - # Use lowest component temperature_crit as starting point - # Starting high and moving down generally works better, - # as it under-predicts next step due to exponential form of - # Psat. - # Subtract 1 to avoid potential singularities at Tcrit - Tbub_initial = ( - min(blk.params.get_component(j).temperature_crit.value for j in raoult_comps) - - 1 - ) - - m = Block() - blk.add_component("tbub_initialise", m) - m.Tbub = Var(initialize=Tbub_initial) - m.f = Expression(expr=( - sum( - get_method(blk, "pressure_sat_comp", j)( - blk, blk.params.get_component(j), m.Tbub * T_units - ) - * blk.mole_frac_comp[j] - for j in raoult_comps - ) - + sum( - blk.mole_frac_comp[j] - * blk.params.get_component(j) - .config.henry_component[liquid_phase]["method"] - .return_expression(blk, liquid_phase, j, m.Tbub * T_units) - for j in henry_comps - ) - - blk.pressure - )) - m.df = Expression(expr=( - sum( - get_method(blk, "pressure_sat_comp", j)( - blk, blk.params.get_component(j), m.Tbub * T_units, dT=True - ) - * blk.mole_frac_comp[j] - for j in raoult_comps - ) - + sum( - blk.mole_frac_comp[j] - * blk.params.get_component(j) - .config.henry_component[liquid_phase]["method"] - .dT_expression(blk, liquid_phase, j, m.Tbub * T_units) - for j in henry_comps - ) - )) - m.tbub_constraint = Constraint(expr=( - m.Tbub == m.Tbub - m.f / m.df - )) - solver = SolverFactory("ipopt") - solver.options["tol"] = 1e-8 - solver.solve(m) - - result = m.Tbub.value - blk.del_component(m) # cleanup - return result - - -def _custom_estimate_Tdew(blk, T_units, raoult_comps, henry_comps, liquid_phase): - """ - Function to estimate dew point temperature - - Args: - blk: StateBlock to use - T_units: units of temperature - raoult_comps: list of components that follow Raoult's Law - henry_comps: list of components that follow Henry's Law - liquid_phase: name of liquid phase - - Returns: - Estimated dew point temperature as a float. - - """ - # Use lowest component critical temperature - # as starting point - # Subtract 1 to avoid potential singularities at Tcrit - Tdew_initial = ( - min(blk.params.get_component(j).temperature_crit.value for j in raoult_comps) - - 1 - ) - - m = Block() - blk.add_component("tdew_initialise", m) - m.Tdew = Var(initialize=Tdew_initial) - m.f = Expression(expr=( - blk.pressure - * ( - sum( - blk.mole_frac_comp[j] - / get_method(blk, "pressure_sat_comp", j)( - blk, blk.params.get_component(j), m.Tdew * T_units - ) - for j in raoult_comps - ) - + sum( - blk.mole_frac_comp[j] - / blk.params.get_component(j) - .config.henry_component[liquid_phase]["method"] - .return_expression(blk, liquid_phase, j, m.Tdew * T_units) - for j in henry_comps - ) - ) - - 1 - )) - m.df = Expression(expr=-( - blk.pressure - * ( - sum( - blk.mole_frac_comp[j] - / get_method(blk, "pressure_sat_comp", j)( - blk, blk.params.get_component(j), m.Tdew * T_units - ) - ** 2 - * get_method(blk, "pressure_sat_comp", j)( - blk, - blk.params.get_component(j), - m.Tdew * T_units, - dT=True, - ) - for j in raoult_comps - ) - + sum( - blk.mole_frac_comp[j] - / blk.params.get_component(j) - .config.henry_component[liquid_phase]["method"] - .return_expression(blk, liquid_phase, j, m.Tdew * T_units) - ** 2 - * blk.params.get_component(j) - .config.henry_component[liquid_phase]["method"] - .dT_expression(blk, liquid_phase, j, m,Tdew * T_units) - for j in henry_comps - ) - ) - )) - m.Tdew_constraint = Constraint(expr=( - m.Tdew == m.Tdew - m.f / m.df - )) - solver = SolverFactory("ipopt") - solver.options["tol"] = TOL - solver.solve(m) - - result = m.Tdew.value - blk.del_component(m) # cleanup - return result diff --git a/property_packages/modular/modular_extended.py b/property_packages/modular/modular_extended.py index 6e8670e..29cf7f4 100644 --- a/property_packages/modular/modular_extended.py +++ b/property_packages/modular/modular_extended.py @@ -1,703 +1,92 @@ -# extends the IDAES helmholtz property package to include additional properties and methods. -from idaes.models.properties.general_helmholtz.helmholtz_state import HelmholtzStateBlockData, _StateBlock -from idaes.models.properties.general_helmholtz.helmholtz_functions import HelmholtzParameterBlockData -from idaes.models.properties.modular_properties.base.generic_property import GenericParameterBlock, _GenericStateBlock, GenericParameterData, GenericStateBlockData +from pyomo.environ import Block, Constraint +from pyomo.core.base.expression import ScalarExpression +from pyomo.core.base.var import ScalarVar, _GeneralVarData, VarData from idaes.core import declare_process_block_class -from property_packages.utils.add_extra_expressions import add_extra_expressions -from pyomo.environ import Constraint, Block -from pyomo.core.base.expression import Expression, ScalarExpression, _GeneralExpressionData, ExpressionData -from pyomo.core.base.var import IndexedVar, ScalarVar, Var, _GeneralVarData,VarData -import idaes.logger as idaeslog -from pyomo.environ import SolverFactory -from idaes.core.util.model_statistics import degrees_of_freedom - -from pyomo.environ import ( - Block, - check_optimal_termination, - Constraint, - exp, - Expression, - log, - Set, - Param, - value, - Var, - units as pyunits, - Reference, -) -from pyomo.common.config import ConfigBlock, ConfigDict, ConfigValue, In, Bool -from pyomo.util.calc_var_value import calculate_variable_from_constraint - -# Import IDAES cores -from idaes.core import ( - declare_process_block_class, - PhysicalParameterBlock, - StateBlockData, - StateBlock, - MaterialFlowBasis, - ElectrolytePropertySet, -) -from idaes.core.base.components import Component, __all_components__ -from idaes.core.base.phases import ( - Phase, - AqueousPhase, - LiquidPhase, - VaporPhase, - __all_phases__, -) -from idaes.core.util.initialization import ( - fix_state_vars, - revert_state_vars, - solve_indexed_blocks, -) -from idaes.core.util.model_statistics import ( - degrees_of_freedom, - number_activated_constraints, -) -from idaes.core.util.exceptions import ( - BurntToast, - ConfigurationError, - PropertyPackageError, - PropertyNotSupportedError, - InitializationError, +from idaes.models.properties.modular_properties.base.generic_property import ( + _GenericStateBlock, + GenericParameterData, + GenericStateBlockData, ) -from idaes.core.util.misc import add_object_reference -from idaes.core.solvers import get_solver -import idaes.logger as idaeslog -import idaes.core.util.scaling as iscale -from idaes.core.initialization.initializer_base import InitializerBase - -from idaes.models.properties.modular_properties.base.generic_reaction import ( - equil_rxn_config, -) -from idaes.models.properties.modular_properties.base.utility import ( - get_method, - get_phase_method, - GenericPropertyPackageError, - StateIndex, - identify_VL_component_list, - estimate_Tbub, - estimate_Tdew, - estimate_Pbub, - estimate_Pdew, -) -from idaes.models.properties.modular_properties.phase_equil.bubble_dew import ( - LogBubbleDew, -) -from idaes.models.properties.modular_properties.phase_equil.henry import HenryType +import idaes.models.properties.modular_properties.base.utility as utility +from property_packages.utils.add_extra_expressions import add_extra_expressions -# ------- -from .modular_custom_utils import ( - _init_Tbub, - _init_Tdew, - _init_Pbub, - _init_Pdew -) +# increase max iterations for estimating values for bubble, dew, and +# critical point initialization (ie. temperature_bubble, temperature_dew) +utility.MAX_ITER = 1000 -class _ExtendedGenericStateBlock(_GenericStateBlock): +def fix_state_vars(blk): """ - This class contains methods which should be applied to Property Blocks as a - whole, rather than individual elements of indexed Property Blocks. - """ - def __init__(self, *args, **kwargs): - super().__init__(*args, **kwargs) # Missing argument - - def initialize(blk, *args, **kwargs): - from idaes.models.properties.modular_properties.base import utility - utility.MAX_ITER = 1000 - # return blk._custom_super_initialize(*args, **kwargs) - flags = {} - hold_state = kwargs.get("hold_state", False) - # if not hold_state or True: - # for k, b in blk.items(): - # avaliable_constraints = ["enth_mass","entr_mass","temperature","entr_mol","smooth_temperature","vapor_frac"] - # for name, var in b.define_state_vars().items(): - # fix_var = True - # if name == "flow_mol": - # if hasattr(b.constraints, "flow_mass"): - # fix_var = False - # elif name in ["enth_mol", "pressure"]: - # if not var.is_fixed(): - # # check if any of the constraints exist - # for prop_name in avaliable_constraints: - # # check if it's a constraint that - # if hasattr(b.constraints,prop_name): - # # we don't need to fix the variable - # # but we need to remove this from the list of constraints (it can't be used to define pressure) - # avaliable_constraints.remove(prop_name) - # fix_var = False - # break - # if hasattr(b, prop_name): - # # check if it's a variable that needs to be unfixed - # if type(getattr(b,prop_name)) in (ScalarVar,): - # extra_var = getattr(b,prop_name) - # if extra_var.is_fixed(): - # # we need to move this from the list of constraints, - # # it's used to define the state - # avaliable_constraints.remove(prop_name) - # fix_var = False - # break - # if fix_var: - # for i in var: - # if not var[i].is_fixed(): - # # fix the state variable at its current value - # flags[k, name, i] = True - # var[i].fix() - # from idaes.core.util.initialization import fix_state_vars - # fix_state_vars(blk) - - # blk[0].report_vars("file10.txt") - # blk[0].constraints.deactivate() - # num_vars = 0 - # for c in blk[0].component_data_objects(Var, descend_into=True, active=True): - # num_vars += 1 - # if not c.fixed: - # continue - # print(f"{c}: {c.value}") - # num_constraints = 0 - # for c in blk[0].component_data_objects(Constraint, descend_into=True, active=True): - # num_constraints += 1 - # print(c) - # print("number of constraints", num_constraints) - # print("number of variables", num_vars) - # else: - # blk[0].report_vars("fffile11.txt") + Fix all state variables at their current values, + except if a constraint exists that defines the variable. - # state_vars_fixed = True - kwargs["state_vars_fixed"] = True - # print(degrees_of_freedom(blk[0])) - # for i, v in blk.constraints.item(): - # print(i, v) - # res = - # try: - # blk[0].report_vars("file15.txt") - blk._custom_super_initialize(*args, **kwargs) - # except Exception as e: - # blk[0].report_vars("file13.txt") - # raise Exception("failed?") from e - return flags - # super().release_state(res) # release the state after initialization, so we can just fix the variables we need to fix. + This method is analogous to the fix_state_vars method in + idaes.core.util.initialization, but it allows us to handle + the case where we have constraints that define the state. - # # reactivate the variables that were deactivated - # for i, b in blk.items(): - # for name, value, in deactivated_vars[i].items(): - # getattr(b,name).fix(value) + Args: + blk : An IDAES StateBlock object in which to fix the state variables - # if zero degrees of freedom, resovle the model with these constraints - # Todo: resolve the model with whatever constraints are active, and guess the rest if degrees of freedom > 0 - # for i,b in self.items(): - # if degrees_of_freedom(b) == 0: - # s = SolverFactory('ipopt') - # res = s.solve(b, tee=True) - # print("init:11o",res.solver.termination_condition) - - # for i, b in blk.items(): - # b.constraints.activate() - # flags[i] = {} - # if hold_state: - # # Fix the required state variables for zero degrees of freedom, and return a dictionary of the flags. - # if not hasattr(b.constraints,"flow_mass") and not b.flow_mol.is_fixed(): - # # We need to fix the flow_mol variable - # flags[i]["flow_mol"] = True - # b.flow_mol.fix() - # avaliable_constraints = ["enth_mass","entr_mass","temperature","entr_mol","smooth_temperature","vapor_frac"] - # if not b.enth_mol.is_fixed(): - # # check if any of the constraints exist - # found_constraint = False - # for prop_name in avaliable_constraints: - # # check if it's a constraint that - # if hasattr(b.constraints,prop_name): - # # we don't need to fix the variable - # # but we need to remove this from the list of constraints (it can't be used to define pressure) - # avaliable_constraints.remove(prop_name) - # found_constraint = True - # break - # if hasattr(b,prop_name): - # # check if it's a variable that needs to be unfixed - # if type(getattr(b,prop_name)) in (ScalarVar,): - # var = getattr(b,prop_name) - # if var.is_fixed(): - # # we need to move this from the list of constraints, - # # it's used to define the state - # avaliable_constraints.remove(prop_name) - # found_constraint = True - # break + Returns: + A dict keyed by block index, state variable name (as defined by + define_state_variables) and variable index indicating the fixed status + of each variable before the fix_state_vars method was applied. + """ - # if not found_constraint: - # # we need to fix the variable - # flags[i]["enth_mol"] = True - # print("fixing enthalpy") - # b.enth_mol.fix() - # if not b.pressure.is_fixed(): - # # check if any of the constraints exist - # found_constraint = False - # for prop_name in avaliable_constraints: - # if hasattr(b.constraints,prop_name): - # # we don't need to fix the variable - # # but we need to remove this from the list of constraints (it can't be used to define pressure) - # avaliable_constraints.remove(prop_name) - # found_constraint = True - # break - # if hasattr(b,prop_name): - # # check if it's a variable that needs to be unfixed - # if type(getattr(b,prop_name)) in (ScalarVar,): - # var = getattr(b,prop_name) - # if var.is_fixed(): - # # we need to move this from the list of constraints, - # # it's used to define the state - # avaliable_constraints.remove(prop_name) - # found_constraint = True - # break - # if not found_constraint: - # # we need to fix the variable - # flags[i]["pressure"] = True - # b.pressure.fix() - - # n_cons = 0 - # dof = 0 - # skip = False - # for k in blk.values(): - # if degrees_of_freedom(k) < 0: - # # Skip solve if DoF < 0 - this is probably due to a - # # phase-component flow state with flash - # skip = True - # n_cons += number_activated_constraints(k) - # dof += degrees_of_freedom(k) - # if n_cons > 0 and not skip: - # if dof > 0: - # raise InitializationError( - # f"{blk.name} Unexpected degrees of freedom during " - # f"initialization at property initialization step: {dof}." - # ) - # with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: - # res = solve_indexed_blocks(opt, [blk], tee=slc.tee) - # init_log.info( - # "Property initialization: {}.".format(idaeslog.condition(res)) - # ) - # else: - # raise Exception() - - return flags - - def release_state(self, flags, outlvl=idaeslog.NOTSET): - for i, v in self.items(): - if i not in flags: - continue - for key in flags[i]: - getattr(v,key).unfix() + flags = {} + for k, b in blk.items(): + available_constraints = [ + "enth_mol", + "enth_mass", + "entr_mol", + "entr_mass", + "smooth_temperature", + ] + for name, var in b.define_state_vars().items(): + fix_var = True + if name == "flow_mol": + for prop_name in ["flow_mass", "flow_vol"]: + if hasattr(b.constraints, prop_name): + fix_var = False + break + elif name in ["temperature", "pressure"]: + if not var.is_fixed(): + # check if any of the constraints exist + for prop_name in available_constraints: + if hasattr(b.constraints, prop_name): + # don't fix this variable - it is defined by a constraint + available_constraints.remove(prop_name) + fix_var = False + break + for i in var: + flags[k, name, i] = var[i].is_fixed() + if fix_var: + var[i].fix() + return flags - def _custom_super_initialize( - blk, - state_args=None, - state_vars_fixed=False, - hold_state=False, - outlvl=idaeslog.NOTSET, - solver=None, - optarg=None, - ): - """ - Initialization routine for property package. - Keyword Arguments: - state_args : a dict of initial values for the state variables - defined by the property package. - outlvl : sets output level of initialization routine - optarg : solver options dictionary object (default=None, use - default solver options) - state_vars_fixed: Flag to denote if state vars have already been - fixed. - - True - states have already been fixed by the - control volume 1D. Control volume 0D - does not fix the state vars, so will - be False if this state block is used - with 0D blocks. - - False - states have not been fixed. The state - block will deal with fixing/unfixing. - solver : str indicating which solver to use during - initialization (default = None, use default solver) - hold_state : flag indicating whether the initialization routine - should unfix any state variables fixed during - initialization (default=False). - - True - states variables are not unfixed, and - a dict of returned containing flags for - which states were fixed during - initialization. - - False - state variables are unfixed after - initialization by calling the - release_state method - Returns: - If hold_states is True, returns a dict containing flags for - which states were fixed during initialization. - """ - init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="properties") - solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="properties") - init_log.info("Starting initialization") +class _ExtendedGenericStateBlock(_GenericStateBlock): - res = None + def initialize(blk, *args, **kwargs): + flag_dict = fix_state_vars(blk) - for k in blk.values(): - # Deactivate the constraints specific for outlet block i.e. - # when defined state is False - if k.config.defined_state is False: - try: - k.sum_mole_frac_out.deactivate() - except AttributeError: - pass + # Set state_vars_fixed to True to avoid fixing state variables + # during the initialize method, since this would overdefine + # the block if we are using constraints + kwargs["state_vars_fixed"] = True - if hasattr(k, "inherent_equilibrium_constraint") and ( - not k.params._electrolyte - or k.params.config.state_components == StateIndex.true - ): - k.inherent_equilibrium_constraint.deactivate() + # Call the base class initialize method + super().initialize(*args, **kwargs) - # Fix state variables if not already fixed - if state_vars_fixed is False: - flag_dict = fix_state_vars(blk, state_args) - # Confirm DoF for sanity - for k in blk.values(): - if k.always_flash: - # If not always flash, DoF is probably less than zero - # We will handle this elsewhere - dof = degrees_of_freedom(k) - if dof != 0: - raise BurntToast( - "Degrees of freedom were not zero [{}] " - "after trying to fix state variables. " - "Something broke in the generic property " - "package code - please inform the IDAES " - "developers.".format(dof) - ) + if kwargs.get("hold_state") is True: + return flag_dict else: - # When state vars are fixed, check that DoF is 0 - for k in blk.values(): - if degrees_of_freedom(k) != 0: - # PYLINT-TODO - # pylint: disable-next=broad-exception-raised - raise Exception( - "State vars fixed but degrees of " - "freedom for state block is not zero " - "during initialization." - ) - - # Create solver - opt = get_solver(solver, optarg) - - # --------------------------------------------------------------------- - # If present, initialize bubble, dew , and critical point calculations - for k in blk.values(): - T_units = k.params.get_metadata().default_units.TEMPERATURE - - # List of bubble and dew point constraints - cons_list = [ - "eq_pressure_dew", - "eq_pressure_bubble", - "eq_temperature_dew", - "eq_temperature_bubble", - "eq_mole_frac_tbub", - "eq_mole_frac_tdew", - "eq_mole_frac_pbub", - "eq_mole_frac_pdew", - "log_mole_frac_tbub_eqn", - "log_mole_frac_tdew_eqn", - "log_mole_frac_pbub_eqn", - "log_mole_frac_pdew_eqn", - "mole_frac_comp_eq", - "log_mole_frac_comp_eqn", - ] - - # Critical point - with k.lock_attribute_creation_context(): - # Only need to look for one, as it is all-or-nothing - if hasattr(k, "pressure_crit"): - # Initialize critical point properties - _initialize_critical_props(k) - # Add critical point constraints to cons_list - cons_list += k.list_critical_property_constraint_names() - - # Bubble temperature initialization - if hasattr(k, "_mole_frac_tbub"): - _init_Tbub(k, T_units) + blk.release_state(flag_dict) - # Dew temperature initialization - if hasattr(k, "_mole_frac_tdew"): - _init_Tdew(k, T_units) - # Bubble pressure initialization - if hasattr(k, "_mole_frac_pbub"): - _init_Pbub(k) - - # Dew pressure initialization - if hasattr(k, "_mole_frac_pdew"): - _init_Pdew(k) - - # Solve bubble, dew, and critical point constraints - for c in k.component_objects(Constraint): - # Deactivate all constraints not associated with bubble, dew, - # or critical points - if c.local_name not in cons_list: - c.deactivate() - - # If StateBlock has active constraints (i.e. has bubble, dew, or critical - # point calculations), solve the block to converge these - n_cons = 0 - dof = 0 - for k in blk.values(): - n_cons += number_activated_constraints(k) - dof += degrees_of_freedom(k) - if n_cons > 0: - if dof > 0: - raise InitializationError( - f"{blk.name} Unexpected degrees of freedom during " - f"initialization at bubble, dew, and critical point step: {dof}." - ) - with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: - res = solve_indexed_blocks(opt, [blk], tee=slc.tee) - init_log.info( - "Bubble, dew, and critical point initialization: {}.".format( - idaeslog.condition(res) - ) - ) - # --------------------------------------------------------------------- - # Calculate _teq if required - # Using iterator k outside of for loop - this should be OK as we just need - # a valid StateBlockData an assume they are all the same. - if k.params.config.phases_in_equilibrium is not None and ( - not k.config.defined_state or k.always_flash - ): - for k in blk.values(): - for pp in k.params._pe_pairs: - k.params.config.phase_equilibrium_state[pp].calculate_teq(k, pp) - - init_log.info("Equilibrium temperature initialization completed.") - - # --------------------------------------------------------------------- - # Initialize flow rates and compositions - for k in blk.values(): - - k.params.config.state_definition.state_initialization(k) - - if k.params._electrolyte: - if k.params.config.state_components == StateIndex.true: - # First calculate initial values for apparent species flows - for p, j in k.params.apparent_phase_component_set: - calculate_variable_from_constraint( - k.flow_mol_phase_comp_apparent[p, j], - k.true_to_appr_species[p, j], - ) - # Need to calculate all flows before doing mole fractions - for p, j in k.params.apparent_phase_component_set: - sum_flow = sum( - k.flow_mol_phase_comp_apparent[p, jj] - for jj in k.params.apparent_species_set - if (p, jj) in k.params.apparent_phase_component_set - ) - if value(sum_flow) == 0: - x = 1 - else: - x = value(k.flow_mol_phase_comp_apparent[p, j] / sum_flow) - lb = k.mole_frac_phase_comp_apparent[p, j].lb - if lb is not None and x <= lb: - k.mole_frac_phase_comp_apparent[p, j].set_value(lb) - else: - k.mole_frac_phase_comp_apparent[p, j].set_value(x) - elif k.params.config.state_components == StateIndex.apparent: - # First calculate initial values for true species flows - for p, j in k.params.true_phase_component_set: - calculate_variable_from_constraint( - k.flow_mol_phase_comp_true[p, j], - k.appr_to_true_species[p, j], - ) - # Need to calculate all flows before doing mole fractions - for p, j in k.params.true_phase_component_set: - sum_flow = sum( - k.flow_mol_phase_comp_true[p, jj] - for jj in k.params.true_species_set - if (p, jj) in k.params.true_phase_component_set - ) - if value(sum_flow) == 0: - x = 1 - else: - x = value(k.flow_mol_phase_comp_true[p, j] / sum_flow) - lb = k.mole_frac_phase_comp_true[p, j].lb - if lb is not None and x <= lb: - k.mole_frac_phase_comp_true[p, j].set_value(lb) - else: - k.mole_frac_phase_comp_true[p, j].set_value(x) - - # If state block has phase equilibrium, use the average of all - # _teq's as an initial guess for T - if ( - k.params.config.phases_in_equilibrium is not None - and isinstance(k.temperature, Var) - and not k.temperature.fixed - ): - k.temperature.value = value( - sum(k._teq[i] for i in k.params._pe_pairs) / len(k.params._pe_pairs) - ) - - if outlvl > 0: # TODO: Update to use logger Enum - init_log.info("State variable initialization completed.") - - # --------------------------------------------------------------------- - n_cons = 0 - dof = 0 - skip = False - Tfix = {} # In enth based state defs, need to also fix T until later - for k, b in blk.items(): - if b.params.config.phase_equilibrium_state is not None and ( - not b.config.defined_state or b.always_flash - ): - if not b.temperature.fixed: - b.temperature.fix() - Tfix[k] = True - for c in b.component_objects(Constraint): - # Activate common constraints - if c.local_name in ( - "total_flow_balance", - "component_flow_balances", - "sum_mole_frac", - "phase_fraction_constraint", - "mole_frac_phase_comp_eq", - "mole_frac_comp_eq", - ): - c.activate() - if c.local_name == "log_mole_frac_phase_comp_eqn": - c.activate() - for p, j in b.params._phase_component_set: - calculate_variable_from_constraint( - b.log_mole_frac_phase_comp[p, j], - b.log_mole_frac_phase_comp_eqn[p, j], - ) - elif c.local_name == "equilibrium_constraint": - # For systems where the state variables fully define the - # phase equilibrium, we cannot activate the equilibrium - # constraint at this stage. - if "flow_mol_phase_comp" not in b.define_state_vars(): - c.activate() - - for pp in b.params._pe_pairs: - # Activate formulation specific constraints - b.params.config.phase_equilibrium_state[ - pp - ].phase_equil_initialization(b, pp) - - n_cons += number_activated_constraints(b) - dof += degrees_of_freedom(b) - if degrees_of_freedom(b) < 0: - # Skip solve if DoF < 0 - this is probably due to a - # phase-component flow state with flash - skip = True - - if n_cons > 0 and not skip: - if dof > 0: - raise InitializationError( - f"{blk.name} Unexpected degrees of freedom during " - f"initialization at phase equilibrium step: {dof}." - ) - with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: - res = solve_indexed_blocks(opt, [blk], tee=slc.tee) - init_log.info( - "Phase equilibrium initialization: {}.".format(idaeslog.condition(res)) - ) - - # --------------------------------------------------------------------- - # Initialize other properties - for k, b in blk.items(): - for c in b.component_objects(Constraint): - # Activate all constraints except flagged do_not_initialize - if c.local_name not in ( - b.params.config.state_definition.do_not_initialize - ): - c.activate() - if k in Tfix: - b.temperature.unfix() - - # Initialize log-form variables - log_form_vars = [ - "act_phase_comp", - "act_phase_comp_apparent", - "act_phase_comp_true", - "conc_mol_phase_comp", - "conc_mol_phase_comp_apparent", - "conc_mol_phase_comp_true", - "mass_frac_phase_comp", - "mass_frac_phase_comp_apparent", - "mass_frac_phase_comp_true", - "molality_phase_comp", - "molality_phase_comp_apparent", - "molality_phase_comp_true", - "mole_frac_comp", # Might have already been initialized - "mole_frac_phase_comp", # Might have already been initialized - "mole_frac_phase_comp_apparent", - "mole_frac_phase_comp_true", - "pressure_phase_comp", - "pressure_phase_comp_apparent", - "pressure_phase_comp_true", - ] - - for prop in log_form_vars: - if b.is_property_constructed("log_" + prop): - comp = getattr(b, prop) - lcomp = getattr(b, "log_" + prop) - for k2, v in lcomp.items(): - c = value(comp[k2]) - if c <= 0: - c = 1e-8 - lc = log(c) - v.set_value(value(lc)) - - # n_cons = 0 - # dof = 0 - # skip = False - # for k in blk.values(): - # if degrees_of_freedom(k) < 0: - # # Skip solve if DoF < 0 - this is probably due to a - # # phase-component flow state with flash - # skip = True - # n_cons += number_activated_constraints(k) - # dof += degrees_of_freedom(k) - # assert degrees_of_freedom(blk) == 0 - # if n_cons > 0 and not skip: - # if dof > 0: - # raise InitializationError( - # f"{blk.name} Unexpected degrees of freedom during " - # f"initialization at property initialization step: {dof}." - # ) - # blk[0].report_vars(f"file_{blk}.txt") - # # blk[0].temperature.value = 500 - # with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: - # res = solve_indexed_blocks(opt, [blk], tee=slc.tee) - # init_log.info( - # "Property initialization: {}.".format(idaeslog.condition(res)) - # ) - - # --------------------------------------------------------------------- - # Return constraints to initial state - for k in blk.values(): - for c in k.component_objects(Constraint): - if c.local_name in (k.params.config.state_definition.do_not_initialize): - print("activating constraint", c.local_name) - c.activate() - - if res is not None and not check_optimal_termination(res): - raise InitializationError( - f"{blk.name} failed to initialize successfully. Please check " - f"the output logs for more information." - ) - - if state_vars_fixed is False: - if hold_state is True: - return flag_dict - else: - blk.release_state(flag_dict) - - init_log.info( - "Property package initialization: {}.".format(idaeslog.condition(res)) - ) - - -@declare_process_block_class("GenericExtendedStateBlock", block_class=_ExtendedGenericStateBlock) +@declare_process_block_class( + "GenericExtendedStateBlock", block_class=_ExtendedGenericStateBlock +) class GenericExtendedStateBlockData(GenericStateBlockData): def build(self, *args): @@ -706,24 +95,23 @@ def build(self, *args): add_extra_expressions(self) # Add a block for constraints, so we can disable or enable them in bulk self.constraints = Block() - - def constrain(self,name:str,value:float): + + def constrain(self, name: str, value: float): # Value must be a float. TODO: Handle unit conversion. - var = getattr(self,name) + var = getattr(self, name) if type(var) == ScalarExpression: self.constraints.add_component(name, Constraint(expr=var == value)) elif type(var) in (ScalarVar, _GeneralVarData, VarData): var.fix(value) - elif type(var) in ( _GeneralExpressionData, ExpressionData) : - # allowed, but we don't need to fix it (eg. mole_frac_comp in helmholtz) - print(f"Variable {self} {name} is an Expression: {type(var)}") - pass else: - raise Exception(f"Variable {self} {name} is not a Var or Expression: {type(var)}") + raise Exception( + f"Variable {self} {name} is not a Var or Expression: {type(var)}" + ) + @declare_process_block_class("GenericExtendedParameterBlock") class GenericExtendedParameterData(GenericParameterData): def build(self): super().build() # set that we should use the extended state block - self._state_block_class = GenericExtendedStateBlock # type: ignore because it'll get created with declare_process_block_class + self._state_block_class = GenericExtendedStateBlock # type: ignore because it'll get created with declare_process_block_class diff --git a/test.py b/test.py index 04c31ed..8d22f54 100644 --- a/test.py +++ b/test.py @@ -11,6 +11,8 @@ from property_packages.build_package import build_package from idaes.core.util.scaling import get_scaling_factor, set_scaling_factor +print('imported everything') + def report_vars(blk, file): with open(file, "w") as f: for c in blk.component_data_objects(Var, descend_into=True, active=True): @@ -47,30 +49,55 @@ def add_set_bubble_to_blk(blk, val): m = ConcreteModel() m.fs = FlowsheetBlock(dynamic=False) -m.fs.properties = build_package("peng-robinson", ["oxygen", "argon", "nitrogen"]) -m.fs.state = m.fs.properties.build_state_block([0], defined_state=True) +m.fs.properties = build_package("peng-robinson", ["oxygen", "nitrogen"]) +# m.fs.state = m.fs.properties.build_state_block([0], defined_state=True) +m.fs.heater = Heater(property_package=m.fs.properties) +sb1 = m.fs.heater.control_volume.properties_in[0] +sb2 = m.fs.heater.control_volume.properties_out[0] + +# sb1.flow_mol.fix(1) +sb1.constrain("flow_mass", 1) +sb1.mole_frac_comp["oxygen"].fix(1/2) +sb1.mole_frac_comp["nitrogen"].fix(1/2) +sb1.pressure.fix(100000) +sb1.constrain("enth_mol", -10) +sb2.constrain("enth_mol", 2400) +# sb2.temperature.fix(400) + +m.fs.heater.initialize() +solver = SolverFactory("ipopt") +solver.options["max_iter"] = 100 +solver.solve(m, tee=True) -sb = m.fs.state[0] -add_report_vars_to_blk(sb) +for c in sb1.component_data_objects(Block): + print(c) -# report_vars(sb, "file21.txt") -# report_constraints(sb, "file22.txt") -# exit() -sb.flow_mol.fix(1) -sb.mole_frac_comp["oxygen"].fix(0.33) -sb.mole_frac_comp["argon"].fix(0.33) -sb.mole_frac_comp["nitrogen"].fix(0.33) -# sb.enth_mol.fix(30000) -# sb.enth_mol.fix(-9) -# sb.enth_mol.value = 100000 -# sb.enth_mol.fix(4000) -# sb.temperature.fix(410) -sb.temperature.fix(404) -# sb._t1_Vap_Liq.value = 280 -# sb._teq["Vap","Liq"].value = 280 -sb.pressure.fix(200000) +m.fs.heater.report() +exit() -m.fs.state.initialize() +# sb = m.fs.state[0] +# add_report_vars_to_blk(sb) + +# # report_vars(sb, "file21.txt") +# # report_constraints(sb, "file22.txt") +# # exit() +# sb.flow_mol.fix(1) +# sb.mole_frac_comp["oxygen"].fix(0.5) +# # sb.mole_frac_comp["argon"].fix(0.33) +# sb.mole_frac_comp["nitrogen"].fix(0.5) +# # sb.enth_mol.fix(30000) +# # sb.enth_mol.fix(-9) +# # sb.enth_mol.value = 100000 +# sb.constrain("enth_mol", -10) +# # sb.enth_mol.fix(-10) +# # sb.temperature.fix(410) +# # sb.temperature.fix(404) +# # sb._t1_Vap_Liq.value = 280 +# # sb._teq["Vap","Liq"].value = 280 +# sb.pressure.fix(100000) + +# m.fs.state.initialize() +# exit() # sb.enth_mol.fix(60000) # except: # sb.report_vars("file17.txt")