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 ]
