In my experience, Poisson's equation must be "grounded" in at least one location. You might try constraining the value of V on just one boundary.
You might also try solving for a ConvectionTerm on P, rather than a DiffusionTerm on V, e.g. ConvectionTerm(coeff=V_variable.faceGrad, var=P_variable) instead of DiffusionTerm(coeff=P_variable, var=V_variable) Generally, though, I hate the drift-diffusion equations and find them challenging to solve. ________________________________________ From: Fangyuan Jiang <fyji...@uw.edu> Sent: Tuesday, August 13, 2019 1:11 PM To: Guyer, Jonathan E. Dr. (Fed); FIPY Subject: Re: Boundary condition problem Thanks Jon for the quick response. I tried to remove the boundary constrain before, but kept showing "RuntimeError: Factor is exactly singular", I couldn't fix that problem by giving different "dt", and don't know what exactly is the problem. Could you please give me some advice? thanks a lot Jiang On Tue, Aug 13, 2019 at 9:53 AM Guyer, Jonathan E. Dr. (Fed) via fipy <fipy@nist.gov<mailto:fipy@nist.gov>> wrote: Jiang - FiPy is no-flux by default, so if you apply no constraints, P should not flow out. By constraining P to be zero on the exteriorFaces, you guarantee that there will be a flux in our out of the external boundary. - Jon ________________________________________ From: fipy-boun...@nist.gov<mailto:fipy-boun...@nist.gov> <fipy-boun...@nist.gov<mailto:fipy-boun...@nist.gov>> on behalf of Fangyuan Jiang <fyji...@uw.edu<mailto:fyji...@uw.edu>> Sent: Tuesday, August 13, 2019 12:25 PM To: FIPY Subject: Boundary condition problem Good morning! Recently, I start to use fipy to solve partial differential equations, but was stuck in the boundary conditions constrain. I want to constrain the variable P(x,t) within the system, but can't stop it from flowing out of the system. I tried all kinds of boundary condition in fipy, but all failed. Below is the two coupled equations, with the code that I write, please correct me if I am wrong, I would be very grateful, since I don't have any friends who has used fipy to discuss with. Therefore, any of your suggestions would be much appreciated. Best regards Jiang Graduate students of UW Equations: ∇2V(x,t)=-a1* P(x, t) dP(x, t)/dt = a2* ∇(∇V(x,t)*P(x,t)) boundary condition : y01=np.zero(nx) y01[30 :70]=1e21 P(x, t=0) = y01 L=1e-6, nx=100, dx= L/nx, a1= 7.23164*e10, a2=1.1e-13 import numpy as np import matplotlib.pyplot as plt from fipy import Variable, FaceVariable, CellVariable, Grid1D, TransientTerm, DiffusionTerm, Viewer, ConvectionTerm L=1e-6 nx=100 dx= L/nx a1=7.23164e-10 a2=1.1e-13 x = np.linspace(0, L, nx) y01 = np.zeros(nx) y01[30:70] = 1e21 mesh = Grid1D(dx=dx, nx=nx) P_variable = CellVariable(mesh=mesh, name='Positive Charge Density', value= y01) V_variable = CellVariable (mesh=mesh, name='potential', value=0.) Eqn1 = DiffusionTerm(coeff=1, var=V_variable) == -a1 * P_variable Eqn2 = TransientTerm(coeff=1, var=P_variable) == a2 * DiffusionTerm(coeff=P_variable, var=V_variable) V_variable.constrain(0., mesh.exteriorFaces) # P_variable.faceGrad.constrain(0., mesh.exteriorFaces) eq = Eqn1 & Eqn2 steps = 1000 time_step = 0.1 for step in range(steps): eq.solve(dt=time_step) if step%1000==0: viewer = Viewer(vars=(P_variable,), datamin=0, datamax=1.5e21) # plt.pause(1) viewer.plot() _______________________________________________ fipy mailing list fipy@nist.gov<mailto: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 ]
