Hi, Kyle. I'll go ahead and keep this on the list for posterity (and also to allow other people to correct me if I say incorrect things).
On Sun, Nov 23, 2014 at 12:29 AM, Kyle Briton Lawlor <[email protected]> wrote: > Hi, Ray. > > I’ve been taking a look at your code to learn some stuff. I have a few > questions. They are not particularly related to the questions you are > asking, so I thought I would message you separate from the FiPy list. > > Anyways, I am working on a problem where I need to invoke a boundary > condition involving a gradient and I saw in this code you had a similar > thing. > > phi.constrain(0., mesh.facesLeft) > phi.constrain(V, mesh.facesRight) > c_m.constrain(c0, mesh.facesLeft) > c_m.faceGrad.constrain(c_m.faceValue*phi.faceGrad, mesh.facesRight) > > My question is does the value in parentheses here, > (c_m.faceValue)*phi.faceGrad, change as the solution is computed for c_m? I > suppose it is fixed here, c_m.faceValue=1. Is this correct? > The value of c_m.faceValue does indeed change as the solution is computed for c_m. So, no c_m.faceValue is not fixed, nor is phi.faceGrad. Both are terms which are dependent on the governing equations as well as the other boundary conditions, so that the computed answer is the solution such that the normal gradient of c_m at the boundary is equal to the value of c_m at the boundary times the normal gradient of phi at the boundary. You can verify this by post-computing n*grad(c_m), c_m,and n*grad(phi) at the boundary. Then, you can see that n*grad(c_m) = c_m * n*grad(phi), where all of those terms come directly from the solution variables. (here, n* is the external surface normal vector dot operation). This may be made a bit more clear in the way the analytical solution is worked out in the original attachment, and agreement between analytical and numerical results confirms that the non-linear, coupled boundary condition is treated properly in FiPy. This is one of the very nice aspects of FiPy in my opinion. Although convergence is perpetually an issue with coupled non-linear problems, the ease with which you can implement even particularly nasty-looking boundary conditions is extremely convenient. > > Best, > Kyle > > > On Nov 19, 2014, at 5:12 PM, Raymond Smith <[email protected]> wrote: > > Hi, FiPy. > > I'm trying to solve a number of problems involving (quasi-)neutral > electrolyte solutions. I've attached some notes I made a while ago on one > of the simplest example problems, for which there is an analytical > solution. The equations are non-linear, and my typical boundary conditions > are also often non-linear, as in the attached example. > > I'm curious if there are suggestions for improving convergence for this > type of problem in FiPy. In the simple example (code here > <https://gist.github.com/raybsmith/6e724e6f1ecc98f9303b> and below), we > have the analytical solution for any applied voltage. For moderate positive > applied voltages, I can get convergence -- up to around 10 thermal volts. > However, for negative applied voltages (which corresponds to depletion of > the salt in the electrolyte), I only get convergence down to about -0.5 > thermal volts, which is notably less than I'd like to simulate. I realize > that sometimes solving non-linear problems numerically is simply difficult, > but perhaps there are some things I could try that someone here could > recommend? > > I've thought about doing the transient version of the problem and stepping > toward a steady state with more extreme applied potentials. However, > because the electrostatic part is static, one of the equations will always > be stationary, even in the dynamic case. (avoiding electrodynamics here) So > in some similar problems, I've still found I can't apply large potentials. > > Thanks, > Ray > > import fipy as fp > > Lx = 1. > nx = 1000 > dx = Lx/nx > mesh = fp.Grid1D(nx=nx, dx=dx) > phi = fp.CellVariable(name="elec. potential", mesh=mesh, value=0.) > c_m = fp.CellVariable(name="conc.", mesh=mesh, value=1.0) > > #V = -6e-1 > V = -4e-1 > #V = 4.0e0 > c0 = 1. > c_m.setValue(c0) > D_m = 2.03e-9 # m^2/s > D_p = 1.33e-9 # m^2/s > D_p = D_p/D_m > D_m = 1. > > eq1 = (0 == fp.DiffusionTerm(coeff=-D_m, var=c_m) + > fp.DiffusionTerm(coeff=D_m*c_m, var=phi)) > eq2 = (0 == fp.DiffusionTerm(coeff=-(D_p - D_m), var=c_m) + > fp.DiffusionTerm(coeff=-(D_p + D_m)*c_m, var=phi)) > > phi.constrain(0., mesh.facesLeft) > phi.constrain(V, mesh.facesRight) > c_m.constrain(c0, mesh.facesLeft) > c_m.faceGrad.constrain(c_m.faceValue*phi.faceGrad, mesh.facesRight) > > Xcc = mesh.cellCenters[0] > J = 1 - fp.numerix.exp(V) > c_m_analyt = fp.CellVariable(name="analyt. conc.", mesh=mesh) > c_m_analyt.setValue(1 - J*Xcc) > phi_analyt = fp.CellVariable(name="analyt. phi", mesh=mesh) > phi_analyt.setValue(fp.numerix.log(1-J*Xcc)) > > eq = eq1 & eq2 > > res = 1e10 > nIt = 0 > while res > 1e-6: > res = eq.sweep() > print "Iteration:", nIt, "Residual:", res > nIt += 1 > > anstol = 1e-3 > print c_m.allclose(c_m_analyt, rtol=anstol, atol=anstol) > print phi.allclose(phi_analyt, rtol=anstol, atol=anstol) > > > <elyte.pdf>_______________________________________________ > 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 ]
