Are all three methods supposed to be equivalent? I'm not sure what exactly is wrong. Everything seems to be working as far as I can tell. The third method has a gradient of -2.0, which is the initial value of BC_Left and BC_Right as expected. gradExteriorValue is never updated in the loop so it won't be the same as method 1 and 2.
On Wed, Apr 27, 2011 at 9:11 PM, Gyeong-Geun LEE <[email protected]> wrote: > Thanks for your help! > Three methods are implemented. > #1: explicit method > #2: source term + setvalue in loop > #3: source term (revised) > ------------------------------------------------------------ > #!/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 #1", #1 > mesh=mesh, > value=valueInitial) > phi2 = CellVariable(name="source term #2", #2 > mesh=mesh, > value=valueInitial) > phi3 = CellVariable(name="source term revised #3", #3 > mesh=mesh, > value=valueInitial) > # Viewer initiation > viewer = Viewer(vars=(phi,phi2, phi3), datamin=-5., datamax=5.) > # BC for #1 > valueLeft = 1.0 > valueRight = -1.0 > BCs = (FixedValue(faces=mesh.getFacesLeft(), value=valueLeft), > FixedValue(faces=mesh.getFacesRight(), value=valueRight)) > # BC for source term (#2) > BC_Left = Variable() > BC_Left = (phi2[0] - valueLeft) * 2 / dx > BC_Right = Variable() > BC_Right = (valueRight - phi2[-1]) * 2 / dx > grad2 = phi2.getFaceGrad() > flux2 = -D * grad2 > # BCs for source term revised (#3) > grad3 = phi3.getFaceGrad() > gradExteriorValue = FaceVariable(mesh=mesh, rank=1) > gradExteriorValue.setValue(BC_Left, where=mesh.getFacesLeft()) > gradExteriorValue.setValue(BC_Right, where=mesh.getFacesRight()) > Dexterior = FaceVariable(mesh=mesh) > Dexterior.setValue(D, where=mesh.getExteriorFaces()) > Dinterior = FaceVariable(mesh=mesh, value=D) > Dinterior.setValue(0, where=mesh.getExteriorFaces()) > flux3 = -Dinterior * grad3 - Dexterior * gradExteriorValue > # Govern equation > eqX = TransientTerm() == ExplicitDiffusionTerm(coeff=D) # #1 > eqX2 = TransientTerm() == -flux2.getDivergence() # #2 > eqX3 = TransientTerm() == -flux3.getDivergence() # #3 > # 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) > # for a new boundary condition for source term (#2) > grad2setValue(BC_Left, where=mesh.getFacesLeft()) > grad2.setValue(BC_Right, where=mesh.getFacesRight()) > # solve the equations > eqX.solve(var=phi, boundaryConditions=BCs, dt=dt) > eqX2.solve(var=phi2, dt=dt) > eqX3.solve(var=phi3, dt=dt) > if __name__ == '__main__': > viewer.plot() > > > On Wed, Apr 27, 2011 at 10:43 PM, Daniel Wheeler <[email protected]> > wrote: >> >> On Wed, Apr 27, 2011 at 3:50 AM, Gyeong-Geun LEE <[email protected]> wrote: >> > Thanks for your suggestion. >> > However, it did not work. >> > The values at boundary were diverged! >> >> Post your script and I'll take a look. >> >> -- >> Daniel Wheeler >> > > -- Daniel Wheeler
