# 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 <fipy-boun...@nist.gov> on behalf of Fangyuan
> Jiang <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
> 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 ]
```