Thank you, Jonathan I have used your method and it works. After some test, i have another question: in the term "fluxright", is the faceValue implicit or explicit? I doubt it's explicit, because i have adjust some parameters, it needs really small time steps. Is there some suggestion for this?
Thank you very much
------------------------------------------
The adjusted problem is :
Set left side (x=0) concentration to constant, and right side (x=L) is the
"J=Kr*c^2" flux boundary condition.
Then for diffusion equation:
\partial c / \partial t = D\nabla^2 c
c(x,0) = c_0
c(0,t)=c_0
J(L) =K_r*c^2
the steady state analytical solution is :
c(L)=c_1=\frac{\sqrt{b^2+4bc_0}-b}{2}
b=\frac{D}{K_r L}
c(x)=c_0+\frac{c_1-c_0}{L}x
----------------------------------
python code:
from fipy import *
# Mesh
Lx = 1. # Thickness
nx = 100.
mesh = Grid1D(Lx=Lx,nx=nx) # 1D Meshing
c0=1000.0 # initial value
# CellVariables
c = CellVariable(name="Concentration", mesh=mesh,value=c0)
cA = CellVariable(name="Analytical solution", mesh=mesh)
x = mesh.cellCenters[0]
time = Variable() # Time (s)
D = 1.0 # Diffusion Coeffcient (m2/s)
Kr = 1.0 # Recombination rate (m4/s)
# Boundary conditions
c.constrain(c0, mesh.facesLeft) # left side, keep c=c0
fluxRight = (-Kr * c.faceValue**2 * mesh._orientedFaceNormals *
mesh.facesRight).divergence # right side, J=Kr*c**2
# Get Analytical solution
b = D/Kr/Lx
c1 = 0.5*(-b+numerix.sqrt(b**2+4*b*c0)) # Analytical solution of c(x=Lx)
cA.setValue((c1-c0)/Lx*x+c0)
# concentration plots
if __name__ == '__main__':
viewer = Viewer(vars=(c,cA), datamin=0., datamax=1.2*c0)
viewer.plot()
# Diffusion Equation
eqn = TransientTerm(var=c) == DiffusionTerm(coeff=D, var=c) + fluxRight
# Time steps
steps = 500
dt = 2.e-5
while time() < steps*dt:
time.setValue(time() + dt)
eqn.solve(dt=dt, var=c)
print dt,time
if __name__ == '__main__':
viewer.plot() # update the concentration plots
if __name__ == '__main__':
raw_input("Calculation finished. Press <return> to proceed...")
On Fri, Jul 19, 2013 at 3:06 AM, 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 ]
>
--
Long.ZHANG
Fusion Reactor & Materials Division
Southwestern Institute of Physics
Address: P.O.Box 432, Chengdu, Sichuan, 610041, China
Tel: +86-28-82850436; Fax: +86-28-82850956
Email: [email protected] or [email protected]
1D_diffusion_Kr_SS.py
Description: Binary data
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
