Re: Boundary condition problem
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.134901960784314 #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/1 #m^2/V*s > > mu_p =1.1e-09/1 #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(var
Re: Boundary condition problem
umpy 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.134901960784314 #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/1 #m^2/V*s > mu_p =1.1e-09/1 #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) > 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 > 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 > mailto:fipy@nist.gov>> wro
Re: Boundary condition problem
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 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 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> mailto:fipy-boun...@nist.gov>> on behalf of Fangyuan Jiang 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 ]
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> 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 on behalf of Fangyuan > Jiang > 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 ]
Re: Boundary condition problem
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 on behalf of Fangyuan Jiang 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 ]
