### 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 *
> 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 *
> > # Continuity Equation
> > Pion.equation = TransientTerm(coeff=1., var=Pion) == mu_p *
> 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 *
> > # Continuity Equation
> > Nion.equation = TransientTerm(coeff=1., var=Nion) == -mu_n *
> 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)
> # Continuity Equation
> Pion.equation = TransientTerm(coeff=1., var=Pion) == mu_p *
> 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)
> # Continuity Equation
> Nion.equation = TransientTerm(coeff=1., var=Nion) == -mu_n *
> 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.
>   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.
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

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

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)

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

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
>
> 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)
>
> 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

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)

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 ]

```