Units not consistent when differentiate returns zero about pyomo HOT 8 OPEN

Interesting. I've never even considered what happens to units when using differentiate.

@dallan-keylogic - can you provide a little more context about what you are trying to do? Where are these derivatives getting used?

from pyomo.

@michaelbynum , they're getting used as part of building up an expression for excess molar enthalpy in an activity coefficient model. Initially, I was using AD to build up the entire expression, but it generally builds expressions that are considerably larger than those I can generate by hand. Right now, I'm using it to get derivatives of other user-defined properties to avoid requiring the user to provide the derivative in addition to the property.

v = pyunits.convert(
    getattr(b, pname + "_vol_mol_solvent"), pyunits.m**3 / pyunits.mol
)
dv_dT = differentiate(
    v,
    b.temperature,
    mode=Modes.reverse_symbolic
)

eps = getattr(b, pname + "_relative_permittivity_solvent")
d_eps_dT = getattr(b, pname + "_d_relative_permittivity_solvent_dT")
A_DH = getattr(b, pname + "_A_DH")
Ix = getattr(b, pname + "_ionic_strength")
rho = ClosestApproach
# Temperature derivative of Debye Huckel term
dA_dT = -A_DH / 2 * (
    1/v * dv_dT + 3 * (1/eps) * d_eps_dT
)

from pyomo.

Could you use numeric differentiation instead of symbolic? Or you need symbolic?

from pyomo.

This is in the context of building up a big Pyomo expression, so it has to be symbolic.

from pyomo.

I think we could add some logic to the differentiate function that checks if an int or a float is being returned and, if so, makes it a Pyomo expression with units that match what is expected. Thoughts, @jsiirola?

from pyomo.

Oh, this is actually because of the way we initialize the derivative map here:

from pyomo.

No, I think this is probably easier to handle inside of differentiate.

from pyomo.

Somewhat tangential to this issue, I have another entry in the "zero should be treated as being equivalent with any set of units" file:

pyunits.convert(
      N_A * (q_e)**2/(8 * pi * eps_vacuum * eps_solvent)
      * (1 + b.temperature / eps_solvent * d_eps_solvent_dT)
      * sum(
          b.mole_frac_phase_comp_true[pname, s]
          * abs(cobj(b, s).config.charge) ** 2
          / cobj(b, s).born_radius
          for s in b.params.ion_set
      ),
      to_units=units.ENERGY/units.AMOUNT
)

I have this expression as part of an eNRTL expression, and have just debugged an issue arising when b.params.ion_set is empty. The summation, instead of picking up units of 1/length, returns the integer zero instead, which leads to dimensional inconsistency when the unit conversion is performed.

There's a straightforward solution in this case, however: return zero with correct units if len(b.params.ion_set) == 0, otherwise return this expression.

from pyomo.