Repair erroneous units for material_holdup in CV0D when mass basis de…

…fined (#1460)

* repair erroneous units for material_holdup when mass basis defined; add tests

* blk

* add holdup_units = None

* condense changes under has_holdup conditional and only add the change needed for mass flow, preserving all previous use cases

* blk

* double-check 1d and cover other conditional for holdup units in 0d CV

* resolve unit inconsistency when other material flow basis specified and test

* fix a couple of typos in flowsheet_model

* scale the new vars on  _physicalparameterblock

* blk

* remove mods related to other mass flow basis

* cleanup

* add 1d test to check dynamic with mass flow basis as well as unit consistency

---------

Co-authored-by: Andrew Lee <[email protected]>

idaes/core/base/control_volume0d.py

@@ -341,13 +341,30 @@ def _add_material_balance_common(
 
         # Material holdup and accumulation
         if has_holdup:
+            if (
+                self.properties_in[
+                    self.flowsheet().time.first()
+                ].get_material_flow_basis()
+                == MaterialFlowBasis.mass
+            ):
+                holdup_units = units("mass")
+            elif (
+                self.properties_in[
+                    self.flowsheet().time.first()
+                ].get_material_flow_basis()
+                == MaterialFlowBasis.molar
+            ):
+                holdup_units = units("amount")
+            else:
+                holdup_units = None
+
             self.material_holdup = Var(
                 self.flowsheet().time,
                 pc_set,
                 domain=Reals,
                 initialize=1.0,
                 doc="Material holdup in control volume",
-                units=units("amount"),
+                units=holdup_units,
             )
         if dynamic:
             self.material_accumulation = DerivativeVar(

idaes/core/base/flowsheet_model.py

@@ -322,10 +322,10 @@ def _setup_dynamics(self):
         elif self.config.time_units is None and self.config.dynamic:
             raise ConfigurationError(
                 f"{self.name} - no units were specified for the time domain. "
-                f"Units must be be specified for dynamic models."
+                f"Units must be specified for dynamic models."
             )
         elif self.config.time_units is None and not self.config.dynamic:
-            _log.debug("No units specified for stady-state time domain.")
+            _log.debug("No units specified for steady-state time domain.")
         elif not isinstance(self.config.time_units, _PyomoUnit):
             raise ConfigurationError(
                 "{} unrecognised value for time_units argument. This must be "

idaes/core/base/tests/test_control_volume_0d.py

@@ -689,6 +689,8 @@ def test_add_phase_component_balances_dynamic():
     assert isinstance(m.fs.cv.material_accumulation, Var)
 
     assert_units_consistent(m)
+    assert_units_equivalent(m.fs.cv.material_holdup, units.mol)
+    assert_units_equivalent(m.fs.cv.material_accumulation, units.mol / units.s)
 
 
 @pytest.mark.unit
@@ -2651,3 +2653,33 @@ def test_reports():
     m.fs.cv.add_momentum_balances(has_pressure_change=True)
 
     m.fs.cv.report()
+
+
+@pytest.mark.unit
+def test_dynamic_mass_basis():
+    m = ConcreteModel()
+    m.fs = Flowsheet(dynamic=True, time_units=units.s)
+    m.fs.pp = PhysicalParameterTestBlock()
+    m.fs.rp = ReactionParameterTestBlock(property_package=m.fs.pp)
+    m.fs.pp.basis_switch = 2
+    m.fs.cv = ControlVolume0DBlock(
+        property_package=m.fs.pp,
+        reaction_package=m.fs.rp,
+        dynamic=True,
+    )
+
+    m.fs.cv.add_geometry()
+    m.fs.cv.add_state_blocks(has_phase_equilibrium=False)
+    m.fs.cv.add_reaction_blocks(has_equilibrium=False)
+
+    mb = m.fs.cv.add_phase_component_balances()
+
+    assert isinstance(mb, Constraint)
+    assert len(mb) == 8
+    assert isinstance(m.fs.cv.phase_fraction, Var)
+    assert isinstance(m.fs.cv.material_holdup, Var)
+    assert isinstance(m.fs.cv.material_accumulation, Var)
+
+    assert_units_consistent(m)
+    assert_units_equivalent(m.fs.cv.material_holdup, units.kg)
+    assert_units_equivalent(m.fs.cv.material_accumulation, units.kg / units.s)

idaes/core/base/tests/test_control_volume_1d.py

@@ -3850,3 +3850,35 @@ def test_estimate_states_backward():
             )
             assert value(state.pressure) == pytest.approx(2.5e5, rel=1e-8)
             assert value(state.temperature) == pytest.approx(345, rel=1e-8)
+
+
+@pytest.mark.unit
+def test_dynamic_mass_flow_basis_unit_consistency():
+    m = ConcreteModel()
+    m.fs = Flowsheet(dynamic=True, time_units=units.s)
+    m.fs.pp = PhysicalParameterTestBlock()
+    m.fs.pp.basis_switch = 2
+    m.fs.rp = ReactionParameterTestBlock(property_package=m.fs.pp)
+
+    m.fs.cv = ControlVolume1DBlock(
+        dynamic=True,
+        has_holdup=True,
+        property_package=m.fs.pp,
+        reaction_package=m.fs.rp,
+        transformation_method="dae.finite_difference",
+        transformation_scheme="BACKWARD",
+        finite_elements=10,
+    )
+
+    m.fs.cv.add_geometry()
+    m.fs.cv.add_state_blocks(has_phase_equilibrium=False)
+    m.fs.cv.add_reaction_blocks(has_equilibrium=False)
+
+    m.fs.cv.add_phase_component_balances()
+    m.fs.cv.test_var = Var(
+        m.fs.cv.flowsheet().time, m.fs.pp.phase_list, m.fs.pp.component_list
+    )
+
+    assert_units_consistent(m)
+    assert_units_equivalent(m.fs.cv.material_holdup, units.kg / units.m)
+    assert_units_equivalent(m.fs.cv.material_accumulation, units.kg / units.s / units.m)

idaes/core/util/testing.py

@@ -197,6 +197,8 @@ def build(self):
         self.set_default_scaling("material_dens_mol", 113)
         self.set_default_scaling("material_flow_mass", 114)
         self.set_default_scaling("material_dens_mass", 115)
+        self.set_default_scaling("material_flow_dimensionless", 116)
+        self.set_default_scaling("material_dens_dimensionless", 117)
 
     @classmethod
     def define_metadata(cls, obj):
@@ -259,6 +261,8 @@ def build(self):
         self.material_dens_mol = Var(initialize=1, units=units.mol / units.m**3)
         self.material_flow_mass = Var(initialize=1, units=units.kg / units.s)
         self.material_dens_mass = Var(initialize=1, units=units.kg / units.m**3)
+        self.material_flow_dimensionless = Var(initialize=1, units=units.dimensionless)
+        self.material_dens_dimensionless = Var(initialize=1, units=units.dimensionless)
         self.pressure = Var(initialize=1e5, units=units.Pa)
         self.temperature = Var(initialize=300, units=units.K)
 
@@ -280,14 +284,18 @@ def build(self):
     def get_material_flow_terms(b, p, j):
         if b.config.parameters.basis_switch == 2:
             return b.material_flow_mass
-        else:
+        elif b.config.parameters.basis_switch == 1:
             return b.material_flow_mol
+        else:
+            return b.material_flow_dimensionless
 
     def get_material_density_terms(b, p, j):
         if b.config.parameters.basis_switch == 2:
             return b.material_dens_mass
-        else:
+        elif b.config.parameters.basis_switch == 1:
             return b.material_dens_mol
+        else:
+            return b.material_dens_dimensionless
 
     def get_enthalpy_flow_terms(b, p):
         return b.enthalpy_flow