Thank you Ray and Dr. Guyer for the different approaches you suggested. I was not aware of the harmonic mean for face values. Thank you for sharing the references as well.
Best Regards, Aniruddha On Wed, Feb 10, 2016 at 5:01 PM, Guyer, Jonathan E. Dr. < [email protected]> wrote: > 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 ] >
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
