search

usnistgov / fipy

FiPy is a Finite Volume PDE solver written in Python
http://pages.nist.gov/fipy/en/latest
Other
504 stars 148 forks source link

Divide by zero bug in AbstractConvectionTerm #914

Open opensourcefan314 opened 1 year ago

opensourcefan314 commented 1 year ago

AbstractConvectionTerm tries hard to make sure the peclet number is never exactly 0 by replacing hard 0s with +/- 1e20.

However this line can fail to eliminate all zeros: https://github.com/usnistgov/fipy/blob/master/fipy/terms/abstractConvectionTerm.py#L139 diffCoeff = diffCoeff - (diffCoeff == 0) * geomCoeff / pecletLarge The problem occurs when an entry of geomCoeff is 0, or is small enough that geomCoeff / pecletLarge underflows. I'm not sure what the best principled fix would be, but a hacky fix is to follow this line with an additional line diffCoeff = diffCoeff - (diffCoeff == 0) / pecletLarge

Here is a gist that triggers the divide by 0: https://gist.github.com/opensourcefan314/c7fa68bf44e069406f617439a9417518

Executing the gist prints this warning:

/fipy/variables/variable.py:1143: RuntimeWarning: invalid value encountered in true_divide return self._BinaryOperatorVariable(lambda a, b: a / b, other) In this case, the divide by zero seems harmless, in the sense that it doesn't propagate nans or infs to the solver output.