I agree with everything Ray said. In electrochemical systems (which is what Aniruddha was working on last I knew), swapping an arithmetic mean for a harmonic mean makes an enormous difference. In semiconductors, there's an even more non-linear interpolation that's widely used, known as Scharfetter-Gummel.
Another approach that's widely used in phase field is the "double obstacle potential". Functionally, the free energy goes to infinity for concentrations less than zero or greater than one. Algorithmically, you just change the concentration if it goes out of bounds. Mass conservation can be an issue here. Another approach that some people take is to write evolution equations for chemical potentials instead of concentrations. This is widely done in semiconductor simulations. One issue to worry about here is whether you have invertability between concentration and chemical potential (OK for ideal solution, but tricky for anything else): Plapp. Unified derivation of phase-field models for alloy solidification from a grand-potential functional. Phys. Rev. E (2011) vol. 84 (3) pp. 031601 Welland et al. Multicomponent phase-field model for extremely large partition coefficients. Phys Rev E (2014) vol. 89 (1) pp. 012409 A. Choudhury and B. Nestler, Phys. Rev. E 85, 021602 (2012). On Feb 10, 2016, at 2:45 PM, Raymond Smith <[email protected]> wrote: > 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 ] _______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
