Using internal constraint feature in fipy

Question

I am a newbie to fipy. I am trying to solve the following set of coupled pdes ( variables are p,n and psi) using fipy :-

qDn∇2n − qun∇.(n∇ψ) = q(n-10**11)/10**(-6), in Domain 1
qDp∇2p + qup∇.(p∇ψ) = -q(p-10**21)/10**(-6), in Domain 1
e∇2ψ = −(p − n- N) in Domain  1
e∇2ψ = 0 in Domain 2, where Domain 1 is the region defined by y<=h and Domain 2 is the region defined by h<y<=h+tox

subjected to the boundary conditions :

p(x,0)= p0, n(x,0)= n0, psi(x,0)= 0, p(x,h)= p0*exp(-psi/Vth), n(x,h)= n0*exp(psi/Vth), psi(x,h+tox)= 5 

I don't want to solve for the variables 'p' and 'n' in domain 2 but I want to solve for ''psi'' in both the domains. Since fipy does not allow explicit domain partition, I am using the internal constraint feature in fipy to impose the boundary conditions at the position given by y=h , since this physical boundary is actually inside the defined mesh i.e. it is not a mesh boundary. Just for keeping up with the generality, I am willing to define the other boundary conditions at y=0 and y=h+tox also using the internal constraint feature though I suppose that mesh.facesBottom and mesh.facesTop will work for the boundaries- y=0 and y=h+tox respectively. I have already framed the pdes in fipy. But I am facing difficulty in defining the boundary conditions using the internal constraint. The code that I have written has been shown below.

! pip install fipy
! pip install pyparse
from fipy import *

L= 10**(-6)
h= 20**(-6)
tox= 0.1*10**(-6)
q=1.6*10**(-19)
un=0.14
up=0.045
Vth=0.026
Dp= up*Vth
Dn=un*Vth
p0= 10**(21)
n0= 10**(11)
e0=8.854*10**(-12)

mesh1= Grid2D(dx= L/100,nx=100,dy=h/200,ny=200) # mesh for domain 1 for solving for p and n
mesh2= Grid2D(dx= L/100,nx=100,dy=tox/10,ny=10)
mesh3= mesh1+(mesh2+[[0],[h]]) # mesh for Domain 1 and Domain 2 for solving for psi

x,y= mesh3.cellCenters

N= 10**21*(y<=h) # N changes from Domain 1 to Domain 2
e= 11.9*e0*(y<=h)+ 3.9*e0*(y>h) # e changes from Domain 1 to Domain 2


p1=CellVariable(name='hole',mesh=mesh3,hasOld=True,value=p0)
n1=CellVariable(name='electron',mesh=mesh3,hasOld=True,value=n0)
psi=CellVariable(name='potential',mesh=mesh3,hasOld=True,value=1)

p= p1*(y<=h) # for domain partition
n= n1*(y<=h) # for domain partition

mask1=((y==0)) # Boundary 1 
mask2=(y==h) # Boundary2
mask3=(y==(h+tox)) # Boundary 3
largevalue= 1e50

eq1=(DiffusionTerm(coeff=q*Dn,var=n)-ConvectionTerm(coeff=q*un*psi.faceGrad,var=n)==ImplicitSourceTerm(coeff=q/10**(-6),var=n)-q*n0/10**(-6))
eq2=(DiffusionTerm(coeff=q*Dp,var=p)+ConvectionTerm(coeff=q*up*psi.faceGrad,var=p)==ImplicitSourceTerm(coeff=-q/10**(-6),var=p)+q*p0/10**(-6))
eq3=(DiffusionTerm(coeff=e,var=psi)==-q*(p-n-N))

eq= eq1 & eq2 & eq3
eq.solve(vars=(n,p,psi))

I am not being able to impose the boundary conditions using the internal constraint feature in fipy. Any help regarding this would be highly appreciated.

Answer 1

It is not clear what you have tried. Per the link I gave in answer to your previous question, an internal constraint is obtained by adding the terms -ImplicitSourceTerm(coeff=mask * largeValue) + mask * largeValue * value to the equation where you want to constrain the solution variable to value at the location of mask.

If you had tried that, you would have encountered a number of other problems with your script above, none of which have anything to do with internal boundary conditions. It is much easier to help if you provide a Minimal, Reproducible Example.

In succession, the issues I encounter are:

TypeError                                 Traceback (most recent call last)
<ipython-input-105-ff64d9005c3b> in <module>()
     41 
     42 eq= eq1 & eq2 & eq3
---> 43 eq.solve(vars=(n,p,psi))

TypeError: solve() got an unexpected keyword argument 'vars'

solve() can take a single var= argument, but only if you didn't specify the solution variable in the equation. With coupled equations, each Term must specify the variable it applies to and you must not provide a var= argument to solve().

---

TypeError                                 Traceback (most recent call last)
<ipython-input-106-6b1bd62471a5> in <module>()
     41 
     42 eq= eq1 & eq2 & eq3
---> 43 eq.solve()
      :
      :
--> 323             b += numerix.reshape(self.constraintB.ravel(), ids.shape).sum(-2).ravel()
    324 
    325         return (var, L, b)

TypeError: ufunc 'add' output (typecode 'O') could not be coerced to provided output parameter (typecode 'd') according to the casting rule ''same_kind''

This is subtle and occurs because 10**(21) doesn't mean what you think it does. This expression takes the integer 10 and raises it to the 21st power, resulting in the long integer 1000000000000000000000L. FiPy cannot solve for integers. All instances of 10**?? need to be replaced with floating point scientific notation, 1e??.

---

ValueError                                Traceback (most recent call last)
<ipython-input-109-c1eb6812a3e3> in <module>()
     41 
     42 eq= eq1 & eq2 & eq3
---> 43 eq.solve()
      :

ValueError: operands could not be broadcast together with shapes (21000,) (2,2) 

This is also rather subtle. Tuple unpacking like x,y= mesh3.cellCenters assigns numpy ndarrays to x and y, but FiPy needs them to be CellVariables (see issue #325). Instead, do

x = mesh3.cellCenters[0]
y = mesh3.cellCenters[1]

---

TypeError                                 Traceback (most recent call last)
<ipython-input-115-894ff5198dc6> in <module>()
     42 
     43 eq= eq1 & eq2 & eq3
---> 44 eq.solve()
      :
      :
.../fipy/solvers/pysparse/pysparseSolver.pyc in _solve(self)
     77             array /= self.var.unit.factor
     78 
---> 79         self.var[:] = array.reshape(self.var.shape)
      :

TypeError: The value of an `_OperatorVariable` cannot be assigned

This is because you are trying to use the expression p= p1*(y<=h) as a solution variable. FiPy can only solve for plain CellVariables. Solve for p1 (I'd just call it p) and put the constraint in the Poisson equation.

---

At this point, I no longer get errors, but there aren't any boundary conditions. The external BCs are easily done with .contrain(). The internal BCs are as described in the link above. One thing you will encounter is that mask2=(y==h) doesn't do anything because there are no cells with a y centered on h. Simplistically, you can do mask2=((y>h - L/100) * (y < h)) to locate the cells in Domain 1 immediately adjacent to h.

---

Putting it all together,

from fipy import *

L= 1e-6
h= 2e-6
tox= 0.1*1e-6
q=1.6*1e-19
un=0.14
up=0.045
Vth=0.026
Dp= up*Vth
Dn=un*Vth
p0= 1e21
n0= 1e11
e0=8.854*1e-12

mesh1= Grid2D(dx= L/100,nx=100,dy=h/200,ny=200) # mesh for domain 1 for solving for p and n
mesh2= Grid2D(dx= L/100,nx=100,dy=tox/10,ny=10)
mesh3= mesh1+(mesh2+[[0],[h]]) # mesh for Domain 1 and Domain 2 for solving for psi

x = mesh3.cellCenters[0]
y = mesh3.cellCenters[1]

N= 1e21*(y<=h) # N changes from Domain 1 to Domain 2
e= 11.9*e0*(y<=h)+ 3.9*e0*(y>h) # e changes from Domain 1 to Domain 2


p=CellVariable(name='hole',mesh=mesh3,hasOld=True,value=p0)
n=CellVariable(name='electron',mesh=mesh3,hasOld=True,value=n0)
psi=CellVariable(name='potential',mesh=mesh3,hasOld=True,value=1.)

p.constrain(p0, where=mesh3.facesBottom)
n.constrain(n0, where=mesh3.facesBottom)
psi.constrain(0., where=mesh3.facesBottom)
psi.constrain(5., where=mesh3.facesTop)

mask1=((y==0)) # Boundary 1 
mask2=((y>h - L/100) * (y < h)) # Boundary2
mask3=(y==(h+tox)) # Boundary 3
largevalue= 1e50

eq1=(DiffusionTerm(coeff=q*Dn,var=n)-ConvectionTerm(coeff=q*un*psi.faceGrad,var=n)
     == ImplicitSourceTerm(coeff=q/1e-6,var=n)-q*n0/1e-6
     - ImplicitSourceTerm(coeff=mask2 * largevalue, var=n) + mask2 * largevalue * n0*numerix.exp(psi/Vth))
eq2=(DiffusionTerm(coeff=q*Dp,var=p)+ConvectionTerm(coeff=q*up*psi.faceGrad,var=p)
     ==ImplicitSourceTerm(coeff=-q/1e-6,var=p)+q*p0/1e-6
     - ImplicitSourceTerm(coeff=mask2 * largevalue, var=p) + mask2 * largevalue * p0*numerix.exp(-psi/Vth))     
eq3=(DiffusionTerm(coeff=e,var=psi)==-q*(p-n-N)*(y <= h))

eq= eq1 & eq2 & eq3
eq.solve()

This will solve, and it will do something like constraining n and p at h, but the solution still won't be right. The semiconductor drift-diffusion equations are absurdly nonlinear and a single direct .solve() has not hope of working. Sweeping is necessary, but because of aggressively nonlinear expressions like n(x,h)= n0*exp(psi/Vth), successive calls to .sweep() will blow up without damping, Newton-Raphson iterations, and probably some attempt to implement Scharfetter-Gummel discretization. Definitely work in 1D until you get things working satisfactorily and solve problems you know the answer for. In the past, I spent a lot of time trying to model photovoltaics with FiPy [e.g. https://doi.org/10.1063/1.3561487] and I've never bothered to write up an example because it's all so messy and unstable.

For further discussion, I strongly urge you to take this to our mailing list or issue tracker. StackOverflow is not well suited to the back-and-forth needed to troubleshoot a problem like this.