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 <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 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 
<fipy@nist.gov<mailto:fipy@nist.gov>> 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 
> <jpoth...@uw.edu<mailto:jpoth...@uw.edu>> 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 
> <fipy@nist.gov<mailto: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 
> > <jpoth...@uw.edu<mailto:jpoth...@uw.edu>> 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<mailto: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 
> > > <jpoth...@uw.edu<mailto: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<mailto: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<mailto:fipy@nist.gov>
> > http://www.ctcms.nist.gov/fipy
> >   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
> > <Fitted Experimental Data.py>
>
>
> _______________________________________________
> fipy mailing list
> fipy@nist.gov<mailto: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<mailto:fipy@nist.gov>
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
--- Begin Message ---
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 <
fipy@nist.gov> 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 <jpoth...@uw.edu> 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 <
> 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 <jpoth...@uw.edu> 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 <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 ]
> > > <Fitted Experimental Data.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 ]
>

--- End Message ---
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Reply via email to