Scott - Eq2 appears to have a typo:
-Eq2 = DiffusionTerm(coeff=z2*C1, var=Phi) + DiffusionTerm(coeff=1.0, var=C2) +Eq2 = DiffusionTerm(coeff=z2*C2, var=Phi) + DiffusionTerm(coeff=1.0, var=C2) Eq3 has (much) better coupling if written like this: -Eq3 = DiffusionTerm(coeff=z3*C3, var=Phi) + (C3.faceGrad).divergence +Eq3 = DiffusionTerm(coeff=z3*C3, var=Phi) - DiffusionTerm(coeff=-z1/z3, var=C1) - DiffusionTerm(coeff=-z2/z3, var=C2) [Note that this way you don't need to evaluate .faceGrad at the boundaries, so your issue #650 doesn't come up. We don't upwind at the boundaries because it's complicated and expensive to do right and it's unusual to care.] With these changes, the solution converges in one sweep, but does not agree with the Analytical Solution you give from West. I don't understand a couple of things though. As posed, what governs \Phi? I get \Phi = 0 and see no reason why it would be otherwise. Typically, I would expect Poisson's equation to be obeyed. Since you assert charge electroneutrality, this reduces to Laplace's equation, but in that case, the slope of \Phi should be constant (unless the permittivity is non-uniform, but you don't have a permittivity). In the Analytical plot, \Phi has definite curvature. My guess is you don't need Eq3 (it'll get taken care of automatically) and you need some form of Poisson's equation. - Jon > On Jun 11, 2019, at 4:48 PM, Scott Calabrese Barton <[email protected]> wrote: > > I’ve had some trouble using fipy for the 1D, steady-state ternary electrolyte > problem, which involves coupled, nonlinear migration and diffusion. > > A text version of the model is listed below. A more complete development is > at > https://gcc01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fbit.ly%2F2wJDXob&data=02%7C01%7Cjonathan.guyer%40nist.gov%7C1687dc7e080545a63f0708d6eeae58e7%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%7C1%7C636958829930551626&sdata=cSBGKKJG0%2Ft%2FMtTZbnTwenKAVsKPHa%2F1C8NgupkRI0E%3D&reserved=0 > > from fipy import * > > # Parameters > nx=50 > C20=4.0/3.0 > z1, z2, z3= 2.0, -2.0, 1.0 > > # Mesh > mesh = Grid1D(nx=nx, dx=1.0/nx) > > # Variables > C1= CellVariable(mesh = mesh, value = 1.0) > C2= CellVariable(mesh = mesh, value = 4.0/3.0) > Phi = CellVariable(mesh = mesh, value = -0.5) > > # Electroneutrality > C3 = -1/z3 * (z1*C1 + z2*C2) > > # Boundary Conditions > C1.constrain(1.0, mesh.facesRight) > C1.constrain(0.0, mesh.facesLeft) > C2.constrain(C20, mesh.facesRight) > Phi.constrain(0, mesh.facesRight) > > # Governing Equations > Eq1 = DiffusionTerm(coeff=z1*C1, var=Phi) + DiffusionTerm(coeff=1.0, var=C1) > Eq2 = DiffusionTerm(coeff=z2*C1, var=Phi) + DiffusionTerm(coeff=1.0, var=C2) > Eq3 = DiffusionTerm(coeff=z3*C3, var=Phi) + (C3.faceGrad).divergence > > Eqns = Eq1 & Eq2 & Eq3 > > # Solution > res = 1e+10 > restol= 1e-3 > > while res > restol: > res = Eqns.sweep() > print(res) > > > However, this system does not converge. I suspect the issue is with boundary > conditions. I’ve tried zeroing out the diffusion coefficients at the > boundary, but that does not seem to make a difference here. I am able to use > this approach for the analogous binary problem: > https://gcc01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fbit.ly%2F2wOzxfK&data=02%7C01%7Cjonathan.guyer%40nist.gov%7C1687dc7e080545a63f0708d6eeae58e7%7C2ab5d82fd8fa4797a93e054655c61dec%7C1%7C1%7C636958829930551626&sdata=CkRvDWJQl%2FjHRTwtxd6y1Krs3VuALc2bAOXKJYsdmEk%3D&reserved=0 > > Any insights or suggestions would be welcome. > > Cheers, > Scott☘ > > > _______________________________________________ > 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 ]
