Well, I think I spoke a bit too soon. By doing the transient version of the problem and by switching to GMRESSolver, I was able to approach steady state for considerably higher applied voltages (about an order of magnitude higher), which is considerably better and getting quite close to full salt depletion on one side. I still have trouble sweeping to convergence on the time steps for large applied voltages, but perhaps by taking smaller and smaller time steps I'll be able to continue to make progress. Tips for solving the steady state directly (apart from simply beginning with great initial guesses, which isn't feasible for me once I move to 2D and have more complicated boundary conditions) are certainly still welcome, but I think this at least gives me a way to move forward.
Best, Ray On Wed, Nov 19, 2014 at 5:12 PM, Raymond Smith <smit...@mit.edu> 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) > > >
_______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
