PhysicalField mishandles compound units

[achennu]

I am looking to use the units of the PhysicalFields (and thereby Variables)

Here is an example that shows the problem:

>>> from fipy import PhysicalField, Variable
>>> from fipy.tools.dimensions.physicalField import _findUnit
>>> k = PhysicalField(1, '(10**-3*mol/l)**-2/h') 
>>> v = Variable(1, unit='mol/l')

>>>  k
PhysicalField(1,'l**2/h/0.001**2/mol**2')

>>> k.inBaseUnits()  # this converts correctly
PhysicalField(0.000277777777778,'m**6/s/mol**2')

>>> t = v**3 * k  # should have units of mol/l/h
((pow(Variable(value=PhysicalField(1,'mol/l')), 3)) * 1 l**2/h/0.001**2/mol**2)

>>> t.inBaseUnits()
PhysicalField(277777.777778,'mol/s/m**3')

>>> t.inUnitsOf(t.unit)
lib/python2.7/site-packages/fipy/tools/dimensions/physicalField.pyc in conversionTupleTo(self, other)
   1767         """
   1768         if not numerix.alltrue(self.powers == other.powers):
-> 1769             raise TypeError, 'Incompatible units'
   1770 
   1771         # let (s1,d1) be the conversion tuple from 'self' to base units

TypeError: Incompatible units

I thought this may be due to the numeric factor in the units (which I've seen being handled in other cases well), so I factored that out to make the units l**2/mol**2/h. This also leads to a problem, suggesting the issue is more in the parsing of compound units, rather than the numeric factor.

>>> k = PhysicalField(1, 'l**2/mol**2/h')
>>> k
PhysicalField(1,'l**2/h/mol**2')

>>> k.inBaseUnits()
PhysicalField(2.77777777778e-10,'m**6/s/mol**2')

>>> v = Variable(1, unit='mol/l')
>>> t  = v**3 * k
>>> t
((pow(Variable(value=PhysicalField(1,'mol/l')), 3)) * 1 l**2/h/mol**2)
>>> t.unit  # <-- this is wrong, should be mol/l/h
<PhysicalUnit l/h/mol>
>>> t.inBaseUnits().unit  # <-- this is correct
<PhysicalUnit mol/s/m**3>

I also tried to make v into a PhysicalField instead of a Variable, but the problem persists.

This seems like a bug to me, though I can't quite see from the source code where the error occurs yet.

Any help would be appreciated!

[tkphd]

Is it an option to simply define your v as a PhysicalField instead of a Variable? The following code seems to work:

>>> from fipy import PhysicalField                                                            
>>> v = PhysicalField('1 mol/l')                                                              

>>> k = PhysicalField('1 l**2/mol**2/h')
>>> print(k)                                                                                  
1.0 l**2/h/mol**2
>>> print(k.inBaseUnits())
2.77777777778e-10 m**6/s/mol**2

>>> t = v**3 * k                                                                              
>>> print(t)
1.0 mol/h/l

>>> print(t.unit)
<PhysicalUnit mol/h/l>
>>> print(t.inBaseUnits().unit)
<PhysicalUnit mol/s/m**3>

edited for inline returns

[guyer]

There's definitely something wrong in the unit handling for binOps. I'll investigate.

[guyer]

While waiting for pull request #535 to be merged, you should be able to work around this issue with

>>> t = v * v * v * k

as the bug only affects the ** power operator

[achennu]

@guyer Great, thanks for the quick fix. Sorry for the lapse, I was away for a break over the new years!

I have another issue which I've seen mentioned in the mailing list, but have not seen a resolution for yet. Should I open another issue for it?

The creation of meshes (1D in my case) does not seem to take distance units as its input in dx. That is Grid1D(dx=PhysicalField(0.1, 'mm'), nx=100) fails. Why is this so? Without the grid knowing its physical dimension, it seems that having source and diffusion terms with units is rendered futile. A possible workaround is to set the grid in SI units (but without a unit), and then convert all coefficients and terms to base units and 'implicitly' solve it in base units. That sort of works (hard to catch errors), but makes me wonder of the utility of having physical units in the library. Some notes of explanation would be appreciated.

In further progress, a simple TransientTerm() == DiffusionTerm() style equation, throws the error arg 1 must be a 1-dimensional double array of appropriate size. when I try to solve or sweep it. Solver is in PySparse, but other errors occur with other solvers. How can I debug this issue?

[guyer]

@achennu Please file a new ticket. Dimensional meshes are supposed to work, but are evidently never tested, which means they don't work.

[achennu]

@guyer Any ETA on when a next release of fipy (i guess 3.1.4) will be done? I'd like to include fipy as a dependency that is pip installable.

[guyer]

@achennu

1. The next release will happen when all of the issues in https://github.com/usnistgov/fipy/milestone/4 are addressed. Some of those are easy. Some are beyond our control. We're working on it.

2. FiPy is pip installable now, so I don't understand the issue.

3. What does this have to do with units?

[achennu]

Thanks for the info @guyer

Yes, I realize that it is pip installable. The pypi version doesn't include the fixes made in this issue, and requires a special case to do a pip git+git_url install for fipy. Installing it that way includes the fixes for the units. I've noticed an issue with pip installing it. Will report it separately.