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 <[email protected]> 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 <[email protected]<mailto:[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]<mailto:[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]<mailto:[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]<mailto:[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]<mailto:[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]<mailto:[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]<mailto:[email protected]> > > > http://www.ctcms.nist.gov/fipy > > > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > > > > > _______________________________________________ > > fipy mailing list > > [email protected]<mailto:[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]<mailto:[email protected]> > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] _______________________________________________ fipy mailing list [email protected]<mailto:[email protected]> 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 < [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 ] >
--- End Message ---
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
