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 ]
