Hello Jonathan, I ran the code using 3 different methods:
1. Without [] and with printing cL.faceGrad[0][0] 2. With [] and without printing cL.faceGrad[0][0] 3. With the mesh._orientedFaceNormals method you proposed and without printing cL.faceGrad[0][0] I just printed and manually checked the values of cL.faceValue[0], Kr/D*cL.faceValue[0]**2 and Flux/D for all 3 cases for all times and found them to be the same in all 3 cases. This doesn't explain why [] works, but at least it shows all 3 ways are equivalent. Best, Edwin On Thu, Jul 18, 2013 at 2:06 PM, Guyer, Jonathan E. Dr. < [email protected]> wrote: > I don't really know what's going on here. The short answer is, this isn't > how I'd do a flux constraint in the first place. Following > http://www.ctcms.nist.gov/fipy/documentation/USAGE.html#applying-fixed-flux-boundary-conditions > , > > fluxRight = (-Kr * cL.faceValue**2 * mesh._orientedFaceNormals * > mesh.facesRight).divergence > fluxLeft = (-Kr * cL.faceValue**2 * mesh._orientedFaceNormals * > mesh.facesLeft).divergence > > # Diffusion Equation > eqn = TransientTerm(var=cL) == DiffusionTerm(coeff=D, var=cL) + > fluxRight + fluxLeft > > seems to give behavior consistent with what happens your way when you > print cL.faceGrad[0][0]. > > > I don't know why wrapping the constraint in [] even appears to work at > all, much less why it's slow (For the record, the solution is not the same > as when you omit [] and print cL.faceGrad[0][0]). The issue it's attempting > to address is that the flux (and cL.faceGrad) are vectors (rank-1), whereas > -Kr/D*cL.faceValue**2 is a scalar (rank-0). There is a conceptual > inconsistency here. Wrapping the scalar in [] is an attempt to make it > rank-1, but it's not right (it works fine for scalar values (e.g. [[3.]]), > but not for scalar fields; Python and NumPy syntax just don't work that > way). > > One approach would be to make Kr a vector quantity, but since it's labeled > as a recombination *rate*, that seems conceptually incorrect. Instead, I > would guess that what you are asserting is that the normal flux has a > scalar value proportional to the square of the surface concentration, i.e., > > $$\vec{J}\cdot\hat{n} = Kr C^2$$ > > Multiplying the scalar field by a vector field (the face normals), as I > did above, is the right way to get the vector quantity you want. > > > On Jul 18, 2013, at 4:09 AM, ZHANG Long <[email protected]> wrote: > > > Dear Edwin > > > > Thank you for your reply, the [ ] works. But the running speed of the > code has been significantly slowed down. Maybe my laptop is kinda old, but > it's quite obvious. > > > > So we need answers from developers to see what's really going on. > > > > Regards > > Long > > > > δΊ 2013/7/18 16:53, Edwin Sze Lun Khoo ει: > >> Hello Zhang Long, > >> > >> I had the same problem with my code a few days ago and the fix for this > is to add [ ] around the faceGrad constraints. Specifically, try doing: > >> cL.faceGrad.constrain([-Kr/D*cL.faceValue**2], mesh.facesRight) > >> cL.faceGrad.constrain([Kr/D*cL.faceValue**2], mesh.facesLeft) > >> > >> I just ran your code again with this modification and with the print > statement commented out and the code works now. I don't really > know why [ ] needs to be added, so we will have to see what the developers > say. > >> > >> Best, > >> Edwin > >> > >> > >> On Thu, Jul 18, 2013 at 3:15 AM, ZHANG Long <[email protected]> > wrote: > >> Hello,everyone > >> > >> I'm new to fipy, and i am learning it using the examples. There is > weird problem happened: > >> > >> Based on the 1-D diffusion example, i changed the boundary conditions. > >> I need the surface flux is propotional to the square of surface > concentration: J=Kr*C^2. > >> So i use fick's law to calculate the diffusion flux on surface: J = > -D*dC/dx. > >> Then i set the boundary conditions to: > cL.faceGrad.constrain(-Kr/D*cL.faceValue**2, mesh.facesRight) > >> > >> The problem is: > >> In line 46, depending on whether i print the "cL.faceGrad[0][0]" or > not, the calculated results are quite different. What happends when the > "cL.faceGrad[0][0]" is printed? > >> > >> At the same time, i also use the change of total inventory in the > volume to calibrate the flux: J=0.5*d(total inventory)/dt. (i use > cL.cellVolumeAverage*Lx to calculate the total inventory). I found that to > get more precise "faceGrad" on the surface,i need very small mesh grids, is > there better ways to calculate the surface flux? > >> > >> Thank you very much > >> > >> --------------------------------------- > >> Code: > >> > >> from fipy import * > >> from numpy import * > >> > >> # Mesh > >> Lx = 0.5e-3 # Thickness > >> dd = list(2**arange(20))+[2**20]*2 > >> dd+=dd[::-1] > >> nx=len(dd) > >> dx=array(dd)*Lx/sum(dd) > >> mesh = Grid1D(dx=dx) # 1D Meshing > >> > >> # CellVariables > >> cI=1.0e5 > >> cL = CellVariable(name="Lattice concentration", mesh=mesh,value=cI) > >> > >> # Other Variables and constants > >> time = Variable() # Time (s) > >> Temp = 700. # Temperature (K) > >> k = 8.6e-5 # Boltzmann constant (eV/K) > >> D = 2.0e-7 * numerix.exp(-0.4/k/Temp) # Diffusion Coeffcient > (m2/s) > >> Kr=3e-25/numerix.sqrt(Temp)*numerix.exp(2.06/k/Temp) # Recombination > rate (m4/s) > >> > >> Cold=float(cL.cellVolumeAverage*Lx) > >> Flux = 0 > >> > >> # Boundary conditions > >> cL.faceGrad.constrain(-Kr/D*cL.faceValue**2, mesh.facesRight) # J = > -D*dcL/dx = Kr*cL**2 > >> cL.faceGrad.constrain(Kr/D*cL.faceValue**2, mesh.facesLeft) # J = > -D*dcL/dx = Kr*cL**2 > >> > >> # Diffusion Equation > >> eqn = TransientTerm(var=cL) == DiffusionTerm(coeff=D, var=cL) > >> > >> # Time steps > >> steps = 500 > >> dt = 500./steps > >> > >> # concentration plots > >> if __name__ == '__main__': > >> viewer = Viewer(vars=(cL), datamin=0., datamax=1.5*cI) > >> viewer.plot() > >> > >> # Solve the equation > >> while time() < steps*dt: > >> time.setValue(time() + dt) > >> res = 1e25 > >> print time , cL.faceValue[0], Kr/D*cL.faceValue[0]**2, Flux/D > ,cL.faceGrad[0][0] > >> while res/cL[nx/2] > 1e-5: > >> res = eqn.sweep(dt=dt, var=cL) > >> if __name__ == '__main__': > >> viewer.plot() # update the concentration plots > >> Cnew=float(cL.cellVolumeAverage*Lx) # Total inventory > >> Flux = -(Cnew-Cold)/dt/2. > >> Cold=Cnew > >> > >> if __name__ == '__main__': > >> raw_input("Calculation finished. Press <return> to proceed...") > >> > >> > >> _______________________________________________ > >> 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 > >> ] > >> > > > > <ATT00001.c> > > _______________________________________________ > 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 ]
