Hi, Aniruddha. Thanks for the clarification. I've also run into issues with concentrations hitting zero in ideal-solution-like simulations in the past.
One thing worth mentioning, if you're imposing a flux that's too large for the physical system to handle, (e.g. above a diffusion limited current), then the system will certainly "explode" at some finite time when the concentration at the surface reaches zero. The phenomenon of Sand's Time in electrochemistry is related to this issue. So, assuming you're not trying to ask the system to do something unphysical, I've _still_ run into issues in the past often associated with my poor numerical choices. The first simple point is to make sure your time steps aren't too large, such that any given cell isn't overstepping in the direction toward zero. Also, one issue I had at one point was that I was (in an electroneutral electrolyte system with electromigration flux terms) using a linear mean for faceValue's instead of a harmonic mean. That can lead to situations in which flux leaves one side of an "empty" cell with no inlets. Anyway, my first thought is that the issue might be avoided with more careful numerics, but perhaps that is not your issue. I've also run into this issue for more recalcitrant systems where I couldn't figure out how to avoid hitting zero with my ideal-solution-like concentrations. What I ended up doing to "cheat" is not an elegant solution, but it kind-of works. I linearly extended the log terms beyond zero with a piecewise function. That is: f(x) = log(x) if x > epsilon f(x) = log(epsilon) + 1/epsilon * (x - epsilon) which should preserve continuity in the function and first derivative but nothing above that. The result is that concentrations in the system can actually go negative, but they really don't like being there (as long as epsilon is small enough), and at least don't go very negative. Again, not a great solution, but it may be something to try. Ray On Wed, Feb 10, 2016 at 2:30 PM, Aniruddha Jana <[email protected]> wrote: > Hello Ray, > > Thank you for the reply. I should have explained my problem better, sorry > for that. > > Yes, I do not initialize the concentration to zero, unlike in the test > script I have shared. But due to a flux at the boundary in my actual > problem, the concentration in some regions can decrease to very low values > close to zero. > > Any thoughts on how to address the log(x>epsilon) problem that gives the > Runtime Warning, that I mentioned in my previous email? Even though it is a > "warning", the solutions to my equations (not in test script) result in a > "nan" because of undesired values computed by the log function, so I cannot > ignore the Runtime Warning. > > Thanks, > Aniruddha > > On Wed, Feb 10, 2016 at 9:13 AM, Raymond Smith <[email protected]> wrote: > >> Hi, Aniruddha. >> >> I certainly don't know the details of your system, but my first thought >> would be that this term in the free energy should (physically, at least) >> prevent the concentration from actually ever reaching zero, as the chemical >> potential diverges when x->0. So I don't understand why you would >> initialize the system with zero concentrations. From my not-super-fresh >> memory of deriving the ideal solution model from a lattice model in Thermo, >> I didn't think it could really apply until you could at least apply >> Stirling's approximation, implying a moderately large number of the species >> in the system. Anyway, is there a strong reason that you need the >> concentrations to be zero? In my experience with such models, when run >> dynamically, they keep concentrations away from zero. >> >> Ray >> >> On Tue, Feb 9, 2016 at 8:11 PM, Aniruddha Jana <[email protected]> wrote: >> >>> Hello Everyone, >>> >>> I am trying to use an ideal chemical free energy formulation for a >>> problem, where I have terms of the form x*log(x). Computing x*log(x) at x = >>> 0 gives the error "RuntimeWarning: divide by zero encountered in log". >>> Even if I try to prevent computation of x*log(x) at x = 0 by doing (for >>> example): >>> >>> epsilon = 1.0e-3 >>> G.setValue((x>epsilon)*log(x>epsilon)) >>> >>> or >>> >>> G.setValue(x*log(x), where=(x>0)) >>> >>> the error persists. >>> >>> A test script is attached with this email. >>> >>> Looking up online, I read that "log" is computed before "where". Another >>> suggestion says to catch and handle the exception (links below). >>> >>> >>> http://stackoverflow.com/questions/25087769/runtimewarning-divide-by-zero-error-how-to-avoid-python-numpy >>> >>> >>> http://stackoverflow.com/questions/13497891/python-getting-around-division-by-zero >>> >>> Any thoughts on how to get rid of the problem, especially for FiPy, >>> and/or the ideal free energy formulation? >>> >>> Thank You! >>> >>> Best regards, >>> Aniruddha >>> >>> >>> >>> _______________________________________________ >>> fipy mailing list >>> [email protected] >>> http://www.ctcms.nist.gov/fipy >>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] >>> >>> >> >> _______________________________________________ >> fipy mailing list >> [email protected] >> http://www.ctcms.nist.gov/fipy >> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] >> >> > > _______________________________________________ > fipy mailing list > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > >
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
