Phi = 0 is a solution. Are you suggesting that this set of equations has more than one solution?
I realize you are not solving Poisson's equation. My question was to consider the broader physics of the problem. Physically, non-uniform electric field means either that charge is non-uniform, which is not in your problem, or permittivity is non-uniform, which is not in your problem. So, what's in your problem that resolves that physical inconsistency?. > On Jun 14, 2019, at 10:53 AM, Scott Calabrese Barton <s...@msu.edu> wrote: > > To answer your questions first, the non-zero gradient in Phi counteracts the > gradients in C2 and C3 such that the migration and diffusion flux cancel and > overall flux is zero for those species. Regarding the curvature, you’re > thinking of the Poisson equation that is only valid in this context when > concentrations are constant. The variations in concentration give rise to > nonlinear variations in Phi. > > Perhaps I should point out that we can get the numerical solution using other > solvers such as solve_bvp from scipy, and an analytical solution is also > available (that’s the source of the plot). So the problem is well defined. > I’ll be happy to share those calculations if needed. I’m using this problem > as an example to learn more about fipy. > >> On Jun 14, 2019, at 9:10 AM, Guyer, Jonathan E. Dr. (Fed) via fipy >> <fipy@nist.gov> wrote: >> >> I had the same initial thought, but you're not adding new information by >> doing that. You can see this by looking at the matrix stencil that results >> from these equations: >> >> X X C1 >> X C2 >> XX C3 = 0 >> XXX Phi >> >> That bottom right sector of the matrix is empty and so the matrix solvers >> cannot find a solution because you're now trying to solve for four unknowns >> with what are really only three equations. >> >> Things to consider: >> >> In the solution you posted from West, Phi does not have zero gradient at >> either end of the domain. Why? What boundary condition is driving that? >> >> Phi is curved, which implies that there is net charge in the system, but you >> assert zero charge throughout. How does West reconcile that inconsistency? >> >>> On Jun 13, 2019, at 11:25 AM, Scott Calabrese Barton <s...@msu.edu> wrote: >>> >>> That’s true, so we can make C3 explicit (full code at >>> http://bit.ly/2F7UQgR= ): >>> >>> # Variables >>> C1= CellVariable(mesh = mesh, value = 1.0) >>> C2= CellVariable(mesh = mesh, value = 4.0/3.0) >>> Phi = CellVariable(mesh = mesh, value = 0) >>> >>> C30 = -1/z3 * (z1*1.0 + z2*C20) # Electroneutrality >>> C3= CellVariable(mesh = mesh, value = C30) >>> >>> # Boundary Conditions >>> C1.constrain(1.0, mesh.facesRight) >>> C2.constrain(C20, mesh.facesRight) >>> C3.constrain(C30, mesh.facesRight) >>> Phi.constrain(-0.2, mesh.facesRight) >>> C1.constrain(0.0, mesh.facesLeft) >>> >>> # Governing Equations >>> Eq1 = + DiffusionTerm(coeff=z1*C1, var=Phi) + DiffusionTerm(coeff=1.0, >>> var=C1) >>> Eq2 = + DiffusionTerm(coeff=z2*C2, var=Phi) + DiffusionTerm(coeff=1.0, >>> var=C2) >>> Eq3 = + DiffusionTerm(coeff=z3*C3, var=Phi) + DiffusionTerm(coeff=1.0, >>> var=C3) >>> >>> Eq4 = ImplicitSourceTerm(coeff=z1, var=C1) + \ >>> ImplicitSourceTerm(coeff=z2, var=C2) + \ >>> ImplicitSourceTerm(coeff=z3, var=C3) >>> >>> Eqns = Eq1 & Eq2 & Eq3 & Eq4 >>> >>> >>> However, this set of equations does not converge. I’m fairly certain that I >>> am mismanaging the boundary conditions. >>> >>> >>>> On Jun 12, 2019, at 5:02 PM, Guyer, Jonathan E. Dr. (Fed) via fipy >>>> <fipy@nist.gov> wrote: >>>> >>>> There's no zero-flux BC on C3 because there's no such thing as C3 as far >>>> as FiPy is concerned. The boundary conditions on C3 are entirely >>>> determined by the boundary conditions on C1 and C2 (Dirichlet on the right >>>> and Robin on the left). The solution obtained satisfies the equations and >>>> boundary conditions given to FiPy. >>>> >>>>> On Jun 12, 2019, at 2:19 PM, Scott Calabrese Barton <s...@msu.edu> wrote: >>>>> >>>>> Jon, thanks for catching that typo. I wish it were the problem! >>>>> >>>>> The solution with grad phi=0 does not satisfy the zero-flux boundary >>>>> condition on species 3. In the absence of a potential gradient, that >>>>> species would need to also have zero gradient (like species 2 in the >>>>> current solution). >>>>> >>>>> So the question becomes how to enforce that zero-flux BC. Any >>>>> suggestions? I could drop the conservation equation and replace it with >>>>> something that enforces zero-flux everywhere, but I’d prefer to use a >>>>> boundary condition as that would be more applicable in higher dimensions. >>>>> >>>>> ☘ >>>>> >>>>> >>>>>> On Jun 12, 2019, at 1:36 PM, Guyer, Jonathan E. Dr. (Fed) via fipy >>>>>> <fipy@nist.gov> wrote: >>>>>> >>>>>> 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 <s...@msu.edu> >>>>>>> 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 http://bit.ly/2wJDXob== >>>>>>> >>>>>>> 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: >>>>>>> http://bit.ly/2wOzxfK== >>>>>>> >>>>>>> Any insights or suggestions would be welcome. >>>>>>> >>>>>>> Cheers, >>>>>>> Scott☘ >>>>>>> >>>>>>> >>>>>>> _______________________________________________ >>>>>>> fipy mailing list >>>>>>> fipy@nist.gov >>>>>>> http://www.ctcms.nist.gov/fipy== >>>>>>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy== >>>>>>> ] >>>>>> >>>>>> >>>>>> _______________________________________________ >>>>>> fipy mailing list >>>>>> fipy@nist.gov >>>>>> http://www.ctcms.nist.gov/fipy== >>>>>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy== ] >>>>>> >>>>> >>>>> >>>>> _______________________________________________ >>>>> fipy mailing list >>>>> fipy@nist.gov >>>>> http://www.ctcms.nist.gov/fipy== >>>>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy== ] >>>> >>>> >>>> _______________________________________________ >>>> fipy mailing list >>>> fipy@nist.gov >>>> http://www.ctcms.nist.gov/fipy== >>>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy== ] >>>> >>> >> >> >> _______________________________________________ >> fipy mailing list >> fipy@nist.gov >> http://www.ctcms.nist.gov/fipy= >> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy= ] >> > _______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
