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: [email protected] <[email protected]> on behalf of Fangyuan Jiang <[email protected]> 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 [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
