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