Dear Jonathan, Thanks for that reply, your suggestions are very helpful. I tried the”sweeping to convergence” method to replace the coupled method to solve the problem, but actually find the result to be almost the same.I think you are right, any tiny changes in the program will lead to enormous change of the whole system, like if I define different dx and nx for the modelling, the results will be totally different. Still thanks for your help, I talked with my professor about your suggestions, we would like to thank you for all your help.
best wishes Jiang On Thu, Aug 15, 2019 at 6:26 AM Guyer, Jonathan E. Dr. (Fed) via fipy < fipy@nist.gov> wrote: > [I'm kicking this back to the list. We don't have the resources to provide > private support.] > > Fangyuan - > > As I discussed last week when Justin Pothoof was asking about the same > issues, your equations are not coupled because Pion.equation only > implicitly depends upon Pion, Nion.equation upon Nion, and > potential.equation upon potential. To be coupled, your equations must > involve terms that have `var=` more than one solution variable. You can > achieve this with ImplicitSourceTerms for the recombination in the drift > diffusion equations and for the charge density in Poisson's equation. > > That said, you do not *have* to solve these equations in a coupled manner; > you can solve them sequentially and sweep them to convergence. > > The equations are not difficult to solve because of the boundary > conditions. They are difficult to solve because tiny changes in Pion and > Nion lead to enormous changes in potential, which, in turn, leads to > enormous electromigration, which causes enormous changes in Pion and Nion. > They're an unstable nightmare. > > Look into Scharfetter-Gummel discretization (FiPy includes a > ScharfetterGummelFaceVariable which approximates this algorithm, but it's > not really right, and it's not documented). > > Some people attempting to solve these equations apply physics-based > preconditioners. FiPy allows for preconditioning, but we've never tried to > build the custom preconditioners that seem to be needed for these equations. > > You may also want to read up on pseudo-Fermi level formulations and > Slotboom variables. Using pseudo-Fermi levels, I've had some luck with the > stabilization algorithm of Brown and Lindsay > doi:10.1016/0038-1101(76)90177-5. > > - Jon > > > > On Aug 15, 2019, at 3:07 AM, Fangyuan Jiang <fyji...@uw.edu> wrote: > > > > Dear Jon, > > > > Sorry to disturb you again for my drift-diffusion modelling. I have made > some progress that I would like to share, also maybe you would help me to > check my equations because I still have some problems. > > first of all, I am able to simulate the drift-diffusion equation for > single-ion system, below is the code and it runs well if I define > divergence.(grad.potential(x,t) * Pion(x,t))as > ConvectionTerm(coeff=V_variable.faceGrad, var=P_variable) instead of > DiffusionTerm(coeff=P_variable, var=V_variable), the defaulted no-flux > boundary condition works for the ConvectionTerm but not the DiffusionTerm, > even though I don't know why. However, if I include the positive ions and > the negative ions together in the modelling , the code runs very weirdly. I > also attached the code below. I realized that it might be because that the > equations is not coupled well, but I have no idea how to couple these > equations while maintaining the defaulted boundary conditions (I mean, > still use the ConvectionTerm in the equation instead of DiffusionTerm). > Last time, you mentioned that drift-diffusion equation is very challenging, > is that because we can't maintain the boundary condition and the equati! > ons-coupling at the same time using Fipy?(which may be related with the > default no-flux boundary condition?) > > > > Thanks in advance! > > Best wishes, > > Fangyuan > > > > import numpy as np > > import matplotlib.pyplot as plt > > from fipy import Variable, FaceVariable, CellVariable, Grid1D, > TransientTerm, DiffusionTerm, Viewer, ConvectionTerm > > > > l = 1e-6 # Length of system in m > > nx = 100 # number of grid points > > dx = l / nx > > x = np.linspace(0, l, nx) > > q = 1.602e-19 # Elementary Charge in C > > epsilon_r = 25 # Relative permittivity of system > > epsilon = epsilon_r * 8.854e-12 # Permittivity of system C/V*m > > > > kb = 1.380649e-23 # J/K > > T = 298.0 # K > > f = kb * T / q # Volts > > mu_n = 1.1e-09 / 10000 # m^2/V*s > > mu_p = 1.1e-09 / 10000 # m^2/V*s > > Dn = f * mu_n # m^2/s > > Dp = f * mu_p # m^2/s > > > > y01 = np.zeros(nx) > > y01[0:10] = 1e21 > > > > mesh = Grid1D(dx=dx, nx=nx) > > > > > > Pion = CellVariable(mesh=mesh, name='Positive ion Charge Density', > value= y01) > > > > potential = CellVariable (mesh=mesh, name='potential', value=0.) > > > > Pion.equation = TransientTerm(coeff=1, var=Pion) == mu_p * > ConvectionTerm(coeff=potential.faceGrad, var=Pion)+ Dp* > DiffusionTerm(coeff=1., var=Pion) > > > > potential.equation = DiffusionTerm(coeff=1, var=potential) == > Pion*q/epsilon > > > > potential.constrain(0., where=mesh.exteriorFaces) > > > > > > eq = Pion.equation & potential.equation > > > > steps = 1000 > > timestep = 0.5 > > > > > > for step in range(steps): > > eq.solve(dt=timestep ) > > if step%1000==0: > > viewer = Viewer(vars=(Pion,), datamin=0., datamax=1e21) > > # viewer = Viewer(vars=(potential,), datamin= -0.03, datamax=0. ) > > plt.pause(1) > > viewer.plot() > > > > > > > > > > > > > > import numpy as np > > import matplotlib.pyplot as plt > > from fipy import Variable, FaceVariable, CellVariable, Grid1D, > ExplicitDiffusionTerm, TransientTerm, DiffusionTerm, Viewer, > ImplicitSourceTerm, ConvectionTerm > > from scipy import special > > > > a1,b1,c1,a2,b2,c2 = [1.07114255e+00, 6.50014631e+05, -4.51527221e+00, > 1.04633414e+00, > > 1.99708312e+05, -1.52479293e+00] > > > > > > q = 1.602e-19 #Elementary Charge > > pini = 154.1581560721245/q #m^-3 > > nini = 134.95618729/q #m^-3 > > > > k1 = 1.8 > > p1 = 17 > > k2 = 17 > > p2 = 1.8 > > > > l = 0.0000134901960784314 #Length of system in m > > nx = 134 #Number of cells in system > > dx = l/nx #Length of each cell in m > > x = np.linspace(0,l,nx) #Array to calculate initial values in > functions below > > > > > > epsilon_r = 25 #Relative permittivity of system > > epsilon = epsilon_r*8.854e-12 #Permittivity of system C/V*m > > kb = 1.38e-23 #J/K > > T = 298 #K > > f = kb*T/q #Volts > > mu_n =1.1e-09/10000 #m^2/V*s > > mu_p =1.1e-09/10000 #m^2/V*s > > Dn = f * mu_n #m^2/s > > Dp = f * mu_p #m^2/s > > > > # k_rec = q*(mu_n+mu_p)/(2*epsilon) > > k_rec = 0. > > > > # y01 = np.zeros(nx) > > # y01[10:20] = 1e21 > > # > > # y02 = np.zeros(nx) > > # y02[110:120] = 1e21 > > > > mesh = Grid1D(dx=dx, nx=nx) #Establish mesh in how many dimensions > necessary > > > > Pion = CellVariable(mesh=mesh, name='Positive ion Charge Density', > value=y01(x)) > > Nion = CellVariable(mesh=mesh, name='Negative ion Charge Density', > value=y02(x)) > > potential = CellVariable(mesh=mesh, name='Potential', value=0.) > > > > #### Equations set-up #### > > > > #In English: dPion/dt = -1/q * divergence.Jp(x,t) - k_rec * Nion(x,t) * > Pion(x,t) where > > # Jp = q * mu_p * E(x,t) * Pion(x,t) - q * Dp * > grad.Pion(x,t) and E(x,t) = -grad.potential(x,t) > > # Continuity Equation > > Pion.equation = TransientTerm(coeff=1., var=Pion) == mu_p * > ConvectionTerm(coeff=potential.faceGrad, > var=Pion)+Dp*DiffusionTerm(coeff=1., var=Pion) > > > > #In English: dNion/dt = 1/q * divergence.Jn(x,t) - k_rec * Nion(x,t) * > Pion(x,t) where > > # Jn = q * mu_n * E(x,t) * Nion(x,t) - q * Dn * > grad.Nion(x,t) and E(x,t) = -grad.potential(x,t) > > # Continuity Equation > > Nion.equation = TransientTerm(coeff=1., var=Nion) == -mu_n * > ConvectionTerm(coeff=potential.faceGrad, > var=Nion)+Dn*DiffusionTerm(coeff=1., var=Nion) > > > > #In English: d^2potential/dx^2 = q/epsilon * Charge_Density and > Charge Density= Pion -Nion > > # Poisson's Equation > > potential.equation = DiffusionTerm(coeff=1., var=potential) == > (q/epsilon)*(Pion-Nion) > > > > > > ### Boundary conditions ### > > # Fipy is defaulted to be no-flux, so we only need to constrain potential > > potential.constrain(0., where=mesh.exteriorFaces) > > > > ### Solve Equations ### > > eq = Pion.equation & Nion.equation & potential.equation > > > > steps = 1000 > > timestep = 0.5 > > for step in range(steps): > > eq.solve(dt=timestep ) > > if (step%1000)==0: > > viewer = Viewer(vars=(potential,), datamin= -2., datamax=2.) > > # viewer = Viewer(vars=(Pion,), datamin=0., datamax=10e21) > > # viewer = Viewer(vars=(Nion,), datamin=0., datamax= 4e21) > > # plt.pause(1) > > viewer.plot() > > # plt.autoscale() > > > > > > On Tue, Aug 13, 2019 at 1:10 PM Guyer, Jonathan E. Dr. (Fed) < > jonathan.gu...@nist.gov> wrote: > > 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 ] >
_______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
