[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 ]

Reply via email to