Thanks for your suggestion. However, it did not work.
The values at boundary were diverged! On Tue, Apr 19, 2011 at 2:07 AM, Daniel Wheeler <[email protected]>wrote: > > On Sun, Apr 17, 2011 at 9:48 PM, Gyeong-Geun LEE <[email protected]> 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 > > >
