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> 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 <fipy-boun...@nist.gov> on behalf of Fangyuan > Jiang <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 > 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 ]
