On Sun, Apr 17, 2011 at 9:48 PM, Gyeong-Geun LEE <gag...@gmail.com> wrote:
> Hi, everyone,
> I am a newbie fipy user. Thanks in advance
> I have a simple 1D grid problem for materials diffusion.
> Because the flux is determined by various non-linear functions, I
> transformed the simple diffusion example to source term for the flux.
> I got the same result between the explicit diffusion and the source term.
> But, I wonder that the boundary condition is right or not.
> Is there a better way of doing that?
I'll try and look through your script and pick out a few things.
> BCs = (FixedValue(faces=mesh.getFacesLeft(), value=valueLeft),
> FixedValue(faces=mesh.getFacesRight(), value=valueRight))
> # BC for source term
> BC_Left = Variable()
> BC_Left = (phi2[0] - valueLeft) * 2 / dx
> BC_Right = Variable()
> BC_Right = (valueRight - phi2[-1]) * 2 / dx
> # Get grad & flux
> grad = phi2.getFaceGrad()
> grad.setValue(BC_Left, where=mesh.getFacesLeft())
> grad.setValue(BC_Right, where=mesh.getFacesRight())
There are better approaches. Does it work? You can split the flux into
exterior and interior parts to clean up the code and avoid having to
update grad in the loop. You don't need to set "BC_Left = Variable()".
That is doing absolutely nothing. You could do this,
>>> gradExteriorValue = FaceVariable(mesh=m, rank=1)
>>> gradExteriorValue.setValue(BC_Left, where=mesh.getFacesLeft())
>>> gradExteriorValue.setValue(BC_Right, where=mesh.getFacesRight())
>>> Dexterior = FaceVariable(mesh=m)
>>> Dexterior.setValue(D, where=mesh.getExteriorFaces())
>>> Dinterior = FaceVariable(mesh=m, value=D)
>>> Dinterior.setValue(0, where=mesh.getExteriorFaces())
>>> flux = -Dinterior * grad - Dexterior * gradExteriorValue
> flux = -D * grad
> # Govern equation
> eqX = TransientTerm() == ExplicitDiffusionTerm(coeff=D)
> eqX2 = TransientTerm() == -flux.getDivergence()
You're trying to solved a coupled system explicitly. In the next
release FiPy will be able to solve this in an implicit manner, which
should bring substantial improvements.
> # Time term
> time = Variable(value=0.0)
> dt = 0.5 * (dx ** 2 / (2 * D))
> print 'dt=', dt
>
> #while time() <= 20:
> for step in range(100):
> time.setValue(time() + dt)
> eqX.solve(var=phi, boundaryConditions=BCs, dt=dt)
> eqX2.solve(var=phi2, dt=dt)
> # for a new boundary condition
> grad.setValue(BC_Left, where=mesh.getFacesLeft())
> grad.setValue(BC_Right, where=mesh.getFacesRight())
> if __name__ == '__main__':
> viewer.plot()
>
> # Final check
> print phi[0:3]
> print phi2[0:3]
> raw_input('Calculation done. Please hit any key...')
>
>
>
--
Daniel Wheeler
