Re: Drift Diffusion Equation Setup
By doing that, you're double counting.
DiffusionTerm(coeff=-mu_p * Pion, var=potental) represents the mathematical
$\nabla\cdot(-\mu_p P_{ion} \nabla V)$
(or divergence.(-mu_p * grad.potential(x,t) * Pion(x,t)) if you prefer).
ConvectionTerm(coeff=-mu_p * potential.faceGrad, var=Pion) *also* represents
the mathematical $\nabla\cdot(-(\mu_p \nabla V) P_{ion}
\equiv \nabla\cdot(-\mu_p P_{ion} \nabla V)$
In DiffusionTerm form, it is implicit in potential. In ConvectionTerm form, it
is implicit in Pion, but either one represents the drift portion of the current
density and both should not be be used together.
From: Justin Pothoof
Sent: Friday, August 9, 2019 8:24 PM
To: Guyer, Jonathan E. Dr. (Fed); FIPY
Subject: Re: Drift Diffusion Equation Setup
When taking the divergence of current density: divergence.(-mu_p *
grad.potential(x,t) * Pion(x,t)) would the divergence not also be applied to
the Pion term, since it is also a function of x. I essentially applied the
product rule to obtain the `Pion * potential.faceGrad.divergence` part, which
in hindsight should have just been a DiffusionTerm. This is how I would derive
the divergence.(-mu_p * grad.potential(x,t) * Pion(x,t)) as =
DiffusionTerm(coeff=-mu_p * Pion, var=potential) + ConvectionTerm(-mu_p *
potential.faceGrad, var = Pion)
Which would just be the drift term for the current density. So, the entire
divergence.(Jp) would be:
DiffusionTerm(coeff=-mu_p * Pion, var=potential) + ConvectionTerm(-mu_p *
potential.faceGrad, var = Pion) - DiffusionTerm(coeff=Dp, var=Pion)
to include the diffusion part.
Please correct me if this is an incorrect understanding, and thank you so much!
Justin Pothoof
The University of Washington - Department of Chemistry
Pre-Candidacy PhD Student
Ginger Group
On Fri, Aug 9, 2019 at 1:36 PM Guyer, Jonathan E. Dr. (Fed) via fipy
mailto:fipy@nist.gov>> wrote:
Justin -
A couple of things:
- Charge Density is not Pion + Nion, it's Pion - Nion
- What are the terms `Pion * potential.faceGrad.divergence` in Pion.equation
and `Nion * potential.faceGrad.divergence` in Nion.equation? I don't believe
they should be there.
- Your equations are really not coupled. Pion.equation is an implicit function
of only Pion, Nion.equation is an implicit function of only Nion, and
potential.equation is an implicit function of only potential. Binding these
equations together with `&` does not gain you anything. Coupling comes from
having implicit parts more than one variable in your equations, e.g., k_rec *
Nion(x,t) * Pion(x,t) could be ImplicitSourceTerm(coeff=k_rec*Nion) in
Pion.equation and ImplicitSourceTerm(coeff=k_rec*Pion) in Nion.equation.
Likewise, divergence.(-mu_p * grad.potential(x,t) * Pion(x,t)) can either be
ConvectionTerm(coeff=-mu_p * potential.faceGrad, var=Pion) or it can be
DiffusionTerm(coeff=-mu_p * Pion, var=potential).
- Boundary conditions in FiPy are no-flux by default, so there's no need to
constrain Pion and Nion.
- Jon
From: Justin Pothoof mailto:jpoth...@uw.edu>>
Sent: Friday, August 9, 2019 2:51 PM
To: Guyer, Jonathan E. Dr. (Fed); FIPY
Subject: Re: Drift Diffusion Equation Setup
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
Re: Drift Diffusion Equation Setup
When taking the divergence of current density: divergence.(-mu_p * grad.potential(x,t) * Pion(x,t)) would the divergence not also be applied to the Pion term, since it is also a function of x. I essentially applied the product rule to obtain the `Pion * potential.faceGrad.divergence` part, which in hindsight should have just been a DiffusionTerm. This is how I would derive the divergence.(-mu_p * grad.potential(x,t) * Pion(x,t)) as = DiffusionTerm(coeff=-mu_p * Pion, var=potential) + ConvectionTerm(-mu_p * potential.faceGrad, var = Pion) Which would just be the drift term for the current density. So, the entire divergence.(Jp) would be: DiffusionTerm(coeff=-mu_p * Pion, var=potential) + ConvectionTerm(-mu_p * potential.faceGrad, var = Pion) *- DiffusionTerm(coeff=Dp, var=Pion) * to include the diffusion part. Please correct me if this is an incorrect understanding, and thank you so much! Justin Pothoof The University of Washington - Department of Chemistry Pre-Candidacy PhD Student Ginger Group On Fri, Aug 9, 2019 at 1:36 PM Guyer, Jonathan E. Dr. (Fed) via fipy < fipy@nist.gov> wrote: > Justin - > > A couple of things: > - Charge Density is not Pion + Nion, it's Pion - Nion > - What are the terms `Pion * potential.faceGrad.divergence` in > Pion.equation and `Nion * potential.faceGrad.divergence` in Nion.equation? > I don't believe they should be there. > - Your equations are really not coupled. Pion.equation is an implicit > function of only Pion, Nion.equation is an implicit function of only Nion, > and potential.equation is an implicit function of only potential. Binding > these equations together with `&` does not gain you anything. Coupling > comes from having implicit parts more than one variable in your equations, > e.g., k_rec * Nion(x,t) * Pion(x,t) could be > ImplicitSourceTerm(coeff=k_rec*Nion) in Pion.equation and > ImplicitSourceTerm(coeff=k_rec*Pion) in Nion.equation. Likewise, > divergence.(-mu_p * grad.potential(x,t) * Pion(x,t)) can either be > ConvectionTerm(coeff=-mu_p * potential.faceGrad, var=Pion) or it can be > DiffusionTerm(coeff=-mu_p * Pion, var=potential). > - Boundary conditions in FiPy are no-flux by default, so there's no need > to constrain Pion and Nion. > > - Jon > > > From: Justin Pothoof > Sent: Friday, August 9, 2019 2:51 PM > To: Guyer, Jonathan E. Dr. (Fed); FIPY > Subject: Re: Drift Diffusion Equation Setup > > 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.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
Re: Drift Diffusion Equation Setup
Justin - A couple of things: - Charge Density is not Pion + Nion, it's Pion - Nion - What are the terms `Pion * potential.faceGrad.divergence` in Pion.equation and `Nion * potential.faceGrad.divergence` in Nion.equation? I don't believe they should be there. - Your equations are really not coupled. Pion.equation is an implicit function of only Pion, Nion.equation is an implicit function of only Nion, and potential.equation is an implicit function of only potential. Binding these equations together with `&` does not gain you anything. Coupling comes from having implicit parts more than one variable in your equations, e.g., k_rec * Nion(x,t) * Pion(x,t) could be ImplicitSourceTerm(coeff=k_rec*Nion) in Pion.equation and ImplicitSourceTerm(coeff=k_rec*Pion) in Nion.equation. Likewise, divergence.(-mu_p * grad.potential(x,t) * Pion(x,t)) can either be ConvectionTerm(coeff=-mu_p * potential.faceGrad, var=Pion) or it can be DiffusionTerm(coeff=-mu_p * Pion, var=potential). - Boundary conditions in FiPy are no-flux by default, so there's no need to constrain Pion and Nion. - Jon From: Justin Pothoof Sent: Friday, August 9, 2019 2:51 PM To: Guyer, Jonathan E. Dr. (Fed); FIPY Subject: Re: Drift Diffusion Equation Setup 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.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)*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 v
Re: Drift Diffusion Equation Setup
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.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)*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:
Re: Drift Diffusion Equation Setup
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 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 > 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 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 > > 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 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
Re: Drift Diffusion Equation Setup
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 < fipy@nist.gov> 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 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 < > 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. > > > > 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 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 > > > Testing.py>___ > > > 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 ] > > > > > ___ > fipy mailing list > fipy@nist.gov >
Re: Drift Diffusion Equation Setup
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 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 > 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 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 > > > Testing.py>___ > > 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 ] > ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
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.
>
> 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 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
> > Testing.py>___
> > 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.subplots_adjust(wspace=0.5,hspace=0.6)
plt.show()___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: Drift Diffusion Equation Setup
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 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 > Testing.py>___ > 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 ]
