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?
GG
-------------------------------------------------------------------
*#!/usr/bin/env python2.5*
*
*
*from fipy import **
*
*
*# mesh*
*nx = 10*
*dx = 1.*
*mesh = Grid1D(nx=nx, dx=dx)*
*
*
*# Diffusivity*
*D = 1.*
*
*
*# Initiation*
*valueInitial = 0.0*
*phi = CellVariable(name="explict diffusion",*
* mesh=mesh,*
* value=valueInitial)*
*phi2 = CellVariable(name="source term",*
* mesh=mesh,*
* value=valueInitial)*
*
*
*# Viewer initiation*
*viewer = Viewer(vars=(phi,phi2), datamin=-5., datamax=5.)*
*
*
*# BC*
*valueLeft = 1.0*
*valueRight = -1.0*
*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())*
*flux = -D * grad*
*
*
*# Govern equation*
*eqX = TransientTerm() == ExplicitDiffusionTerm(coeff=D)*
*eqX2 = TransientTerm() == -flux.getDivergence()*
*
*
*# 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...')*
