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 <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 ] > _______________________________________________ > 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 ]
