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]" i**s 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 ]
