Hi Jon,
I've been continuing to work on my problem. I've implemented the divergence
of the current density mathematically into the positive charge and negative
charge density equations. Now, I'm encountering issues with the boundary
conditions.. the charges should be flowing towards the center of the
system, and ultimately recombine, but the issues I'm seeing is that the
charges are flowing out of the region. I would like to keep all of the
charges confined to the system and prevent them from flowing out of the
left and right boundaries. I've tried implementing neumann boundaries using
Pion.faceGrad.constrain(0., where=mesh.exteriorFaces) and the negative
charge version, but I'm still seeing the charges flow mostly out of the
system. I would appreciate any help!
Thank you,
Justin
import numpy as np
import matplotlib.pyplot as plt
from scipy import special
from fipy import Variable, FaceVariable, CellVariable, Grid1D,
ExplicitDiffusionTerm, TransientTerm, DiffusionTerm, Viewer,
ImplicitSourceTerm, ConvectionTerm
from fipy.tools import numerix
##########################################################
#################''' SET-UP PARAMETERS '''################
##########################################################
a1,b1,c1,a2,b2,c2 = [1.07114255e+00, 6.50014631e+05, -4.51527221e+00,
1.04633414e+00,
1.99708312e+05, -1.52479293e+00]
# Parameters for sum of sines fit (Potential fit)
#a = -3930.03590805
#b, c = 3049.38274411, -4.01434474
# Parameters for exponential fit (Charge Density) Not used yet
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
# Parameters for charge density fit (Susi's fit)
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)*10
#k_rec = 0
#pini*np.exp(a*x)
#nini*np.exp(b*x+c) #Equations for exponential charge
density fits (Not Used Yet)
##################################################################
##############''' INITIAL CONDITION FUNCTIONS '''#################
##################################################################
def y01(x):
"""Initial positive ion charge density"""
return
pini*((special.gamma(k1+p1))/(special.gamma(k1)*special.gamma(p1))*((x/l)**(k1-1))*(1-(x/l))**(p1-1))/7.3572
def y02(x):
""""Initial negative ion charge density"""
return
nini*((special.gamma(k2+p2))/(special.gamma(k2)*special.gamma(p2))*((x/l)**(k2-1))*(1-(x/l))**(p2-1))/7.3572
def y03(x):
"""Initial potential"""
return a1*np.sin(b1*x+c1) + a2*np.sin(b2*x+c2)
mesh = Grid1D(dx=dx, nx=nx) #Establish mesh in how many dimensions necessary
##############################################################################
#################''' SETUP CELLVARIABLES AND EQUATIONS '''####################
##############################################################################
#CellVariable - defines the variables that you want to solve for:
'''Initial value can be established when defining the variable, or
later using 'var.value ='
Value defaults to zero if not defined'''
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=y03(x))
#EQUATION SETUP BASIC DESCRIPTION
'''Equations to solve for each varible must be defined:
-TransientTerm = dvar/dt
-ConvectionTerm = dvar/dx
-DiffusionTerm = d^2var/dx^2
-Source terms can be described as they would appear mathematically
Notes: coeff = terms that are multiplied by the Term.. must be rank-1
FaceVariable for ConvectionTerm
"var" must be defined for each Term if they are not all the
variable being solved for,
otherwise will see "fipy.terms.ExplicitVariableError: Terms
with explicit Variables cannot mix with Terms with implicit
Variables." '''
#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) + Pion *
potential.faceGrad.divergence) + DiffusionTerm(coeff=Dp,var=Pion) -
k_rec*Pion*Nion
#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) + Nion *
potential.faceGrad.divergence) + DiffusionTerm(coeff=Dn,var=Nion) -
k_rec*Pion*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 '''###################
################################################################
Pion.faceGrad.constrain(0., where=mesh.exteriorFaces) #dPion/dx = 0
at the exterior faces of the mesh
Nion.faceGrad.constrain(0., where=mesh.exteriorFaces) #dNion/dx = 0
at the exterior faces of the mesh
potential.constrain(0., where=mesh.exteriorFaces) #potential = 0
at the exterior faces of the mesh
################################################################
#################''' SOLVE EQUATIONS '''########################
################################################################
eq = Pion.equation & Nion.equation & potential.equation #Couple all
of the equations together
steps = 100 #How many time steps to take
dt = 1 #How long each time step is in seconds
if __name__ == "__main__":
#viewer = Viewer(vars=(potential,),datamin=-1.1,datamax=1.1)
#Sets up viewer for the potential with y-axis limits
viewer = Viewer(vars=(Pion,),datamin=0,datamax=1e21) #Sets
up viewer for negative ion density with y-axis limits
#viewer = Viewer(vars=(Nion,),datamin=-1e21,datamax=0) #Sets
up viewer for positive ion density with y-axis limits
for steps in range(steps): #Time loop to step through
eq.solve(dt=dt) #Solves all coupled equation with timestep dt
if __name__ == '__main__':
viewer.plot() #Plots results using matplotlib
plt.pause(1) #Pauses each frame for n amount of time
#plt.autoscale() #Autoscale axes if necessary
Justin Pothoof
The University of Washington - Department of Chemistry
Pre-Candidacy PhD Student
Ginger Group
On Fri, Jul 26, 2019 at 10:50 AM Guyer, Jonathan E. Dr. (Fed) via fipy <
[email protected]> wrote:
> A current density or flux is a rank-1 variable, typically defined on face
> centers. Your expression
>
> J_n = -mu_n * n.harmonicFaceValue * phi.faceGrad + D_n * n.faceGrad
>
> appropriately declares a rank-1 FaceVariable.
>
> If you attempt to assign this value to a CellVariable, FiPy complains
> because face centers don't coincide with cell centers.
>
> Solution variables in FiPy must be CellVariable objects. This is fine,
> because you don't solve for J_n.
>
> The intention of my earlier suggestion was that you should
> *mathematically* combine your expressions for J_n and dn/dt to obtain a 2nd
> order PDE for dn/dt in terms of CellVariables you know, like n, p, and V.
>
> > On Jul 25, 2019, at 5:15 PM, Justin Pothoof <[email protected]> wrote:
> >
> > Great, that makes a lot of sense!
> >
> > I've tried to define the current density as a CellVariable with the
> value J_n = -mu_n * n.harmonicFaceValue * phi.faceGrad + D_n * n.faceGrad
> > as I've seen you describe in the mailing list before. But, I keep
> encountering the error "ValueError: could not broadcast input array from
> shape (135) into shape (134)" with my mesh defined as length of 134.
> >
> > I believe this is caused by the harmonicFaceValue, though I am not sure?
> >
> > Would the following definition for current density also work: J_p.value
> = -mu_p * p * psi.arithmeticFaceValue.divergence + Dn *
> p.aritmeticFaceValue.divergence
> >
> > I apologize for the multiple questions and I'm very grateful for your
> help!
> > Justin
> >
> > On Thu, Jul 25, 2019 at 10:55 AM Guyer, Jonathan E. Dr. (Fed) via fipy <
> [email protected]> wrote:
> > Justin -
> >
> > I would define a function that takes an argument x for each of your
> analytical expressions, e.g.,
> >
> > ```
> > def y01(x):
> > """Initial positive ion charge density"""
> > return
> pini*((special.gamma(k1+p1))/(special.gamma(k1)*special.gamma(p1))*((x/l)**(k1-1))*(1-(x/l))**(p1-1))/7.3572
> >
> > def y02(x):
> > """"Initial negative ion charge density"""
> > return
> nini*((special.gamma(k2+p2))/(special.gamma(k2)*special.gamma(p2))*((x/l)**(k2-1))*(1-(x/l))**(p2-1))/7.3572
> >
> > def y03(x):
> > """Initial potential"""
> > return a1*np.sin(b1*x+c1) + a2*np.sin(b2*x+c2)
> > ```
> >
> > Then you can invoke these functions with either the linspace to generate
> plots like you have, or with the mesh positions, to set the initial
> conditions:
> >
> > ```
> > Pion.value = y01(mesh.x)
> > Nion.value = y02(mesh.x)
> > potential.value = y03(mesh.x)
> > ```
> >
> > - Jon
> >
> > > On Jul 24, 2019, at 1:23 PM, Justin Pothoof <[email protected]> wrote:
> > >
> > > 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 <
> [email protected]> wrote:
> > > Justin -
> > >
> > > What that error means is that if you write 'var=' for any Term, then
> you must write 'var=' for every Term.
> > >
> > > In your equations:
> > >
> > > ```
> > > 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 takes to advise further.
> > >
> > > - Jon
> > >
> > > > On Jul 23, 2019, at 9:09 PM, Justin Pothoof <[email protected]> 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
> > > > [email protected]
> > > > http://www.ctcms.nist.gov/fipy
> > > > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy
> ]
> > >
> > >
> > > _______________________________________________
> > > fipy mailing list
> > > [email protected]
> > > http://www.ctcms.nist.gov/fipy
> > > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
> > > <Fitted Experimental Data.py>
> >
> >
> > _______________________________________________
> > fipy mailing list
> > [email protected]
> > http://www.ctcms.nist.gov/fipy
> > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
> _______________________________________________
> fipy mailing list
> [email protected]
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
