# Re: Drift Diffusion Equation Setup

Thank you Jon,

I will try writing it in one equation as you suggested.  Regarding the
experimental data, I have an initial potential curve described by a sum of
sines fit as well as initial positive/negative charge density curves
described by a specific equation I'll show in a file.

Thanks for the help!
Justin

On Wed, Jul 24, 2019 at 6:06 AM Guyer, Jonathan E. Dr. (Fed) via fipy <
fipy@nist.gov> wrote:

> Justin -
>
> What that error means is that if you write 'var=' for any Term, then you
> must write 'var=' for every Term.
>
>
> 
> Pion.equation = TransientTerm() + k_rec * Pion * Nion +
> ConvectionTerm(coeff=1 / q, var=Jp) == 0
> Nion.equation = TransientTerm() + k_rec * Pion * Nion +
> ConvectionTerm(coeff=-1 / q, var=Jn) == 0
> potential.equation = DiffusionTerm(1 / q * epsilon) + Pion * Nion == 0
> Jp.equation = ImplicitSourceTerm() + ConvectionTerm(coeff=-q * mu_p *
> Pion, var=potential) + ConvectionTerm(coeff=-q * Dp, var=Pion) == 0
> Jn.equation = ImplicitSourceTerm() + ConvectionTerm(coeff=-q * mu_n *
> Nion, var=potential) + ConvectionTerm(coeff=q * Dn, var=Nion) == 0
> 
> FiPy does not know what Variable TransientTerm() applies to in
> Pion.equation and Nion.equation, DiffusionTerm() in potential.equation, nor
> ImplicitSourceTerm() in Jp.equation and Jn.equation.
>
> As a further point, you should not solve Pion.equation and Jp.equation
> separately nor Nion.equation/Jn.equation. Combine them for a single, second
> order PDE each for n and for p. You will want to take care that, e.g., the
> equation should not be taking the gradient (\nabla) of a vector (Jn), which
> would give you a tensor; rather, the expression should be divergence
> (\nabla\cdot) of a vector (Jn), giving a scalar.
>
> As to comparing to your experimental data, I'd need to know what form it
>
> - Jon
>
> > On Jul 23, 2019, at 9:09 PM, Justin Pothoof <jpoth...@uw.edu> wrote:
> >
> > Hello,
> >
> > I understand that modeling the drift diffusion equations are very
> challenging, but I'm having issues actually writing the equations.  I keep
> encountering the error: "fipy.terms.ExplicitVariableError: Terms with
> explicit Variables cannot mix with Terms with implicit Variables."
> >
> > Additionally, I have fitted experimental data that describes what the
> initial conditions for my system should be, but I don't know how to include
> that into FiPy.  I would appreciate any guidance that you can offer.  I
> will include a pdf of what the equations I'm trying to write are as well as
> the file I have written thus far.
> >
> > Thank you,
> > Justin
> > <FiPy
> Testing.py><Equations.pdf>_______________________________________________
> > 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 ]
>

import numpy as np
import matplotlib.pyplot as plt
from scipy import special

a1,b1,c1,a2,b2,c2 = [ 1.04633244e+00,  1.99709309e+03, -1.52480044e+00,
1.07114122e+00,
6.50014924e+03, -4.51527387e+00]
# Parameters for sum of sines fit

pini = 9.634885e14    #Initial peak positive charge density
nini = -8.434762e14   #Initial peak negative charge density
k1 = 1.8
p1 = 17               #Parameters to fit charge density equations
k2 = 17
p2 = 1.8
l = 0.00134901960784314

x = np.linspace(0,l,134)

y01 =
pini*((special.gamma(k1+p1))/(special.gamma(k1)*special.gamma(p1))*((x/l)**(k1-1))*(1-(x/l))**(p1-1))/7.3572
# Initial positive ion charge density

y02 =
nini*((special.gamma(k2+p2))/(special.gamma(k2)*special.gamma(p2))*((x/l)**(k2-1))*(1-(x/l))**(p2-1))/7.3572
# Initial negative ion charge density

y03 = a1*np.sin(b1*x+c1) + a2*np.sin(b2*x+c2)
# Initial potential

plt.figure(1)
plt.subplot(2,1,1)
plt.plot(x,y03)
plt.xlabel('Position / cm')
plt.ylabel('Potential / V')
plt.title('Potential Fit')

plt.subplot(2,2,3)
plt.plot(x,y01)
plt.xlabel('Position / cm')
plt.ylabel('Charge Density / cm^-3')
plt.title('Positive Charge Density Fit',y=1.08)

plt.subplot(2,2,4)
plt.plot(x,y02)
plt.xlabel('Position / cm')
plt.ylabel('Charge Density / cm^-3')
plt.title('Negative Charge Density Fit',y=1.08)
plt.show()
_______________________________________________