Bookmark and Share FiPy: A Finite Volume PDE Solver Using Python
Version 3.1

This Page

Contact

FiPy developers
Jonathan Guyer
Daniel Wheeler
James Warren

Join our mailing list

100 Bureau Drive, M/S 6555
Gaithersburg, MD 20899

301-975-5329 Telephone
301-975-4553 Facsimile

examples.diffusion.mesh20x20ΒΆ

Solve a two-dimensional diffusion problem in a square domain.

This example solves a diffusion problem and demonstrates the use of applying boundary condition patches.

>>> from fipy import *
>>> nx = 20
>>> ny = nx
>>> dx = 1.
>>> dy = dx
>>> L = dx * nx
>>> mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

We create a CellVariable and initialize it to zero:

>>> phi = CellVariable(name = "solution variable",
...                    mesh = mesh,
...                    value = 0.)

and then create a diffusion equation. This is solved by default with an iterative conjugate gradient solver.

>>> D = 1.
>>> eq = TransientTerm() == DiffusionTerm(coeff=D)

We apply Dirichlet boundary conditions

>>> valueTopLeft = 0
>>> valueBottomRight = 1

to the top-left and bottom-right corners. Neumann boundary conditions are automatically applied to the top-right and bottom-left corners.

>>> X, Y = mesh.faceCenters
>>> facesTopLeft = ((mesh.facesLeft & (Y > L / 2))
...                 | (mesh.facesTop & (X < L / 2)))
>>> facesBottomRight = ((mesh.facesRight & (Y < L / 2))
...                     | (mesh.facesBottom & (X > L / 2)))
>>> phi.constrain(valueTopLeft, facesTopLeft)
>>> phi.constrain(valueBottomRight, facesBottomRight)

We create a viewer to see the results

>>> if __name__ == '__main__':
...     viewer = Viewer(vars=phi, datamin=0., datamax=1.)
...     viewer.plot()

and solve the equation by repeatedly looping in time:

>>> timeStepDuration = 10 * 0.9 * dx**2 / (2 * D)
>>> steps = 10
>>> for step in range(steps):
...     eq.solve(var=phi,
...              dt=timeStepDuration)
...     if __name__ == '__main__':
...         viewer.plot()
solution to diffusion problem on a 2D domain with some Dirichlet boundaries

We can test the value of the bottom-right corner cell.

>>> print numerix.allclose(phi(((L,), (0,))), valueBottomRight, atol = 1e-2)
1
>>> if __name__ == '__main__':
...     raw_input("Implicit transient diffusion. Press <return> to proceed...")

We can also solve the steady-state problem directly

>>> DiffusionTerm().solve(var=phi)
>>> if __name__ == '__main__':
...     viewer.plot()
stead-state solution to diffusion problem on a 2D domain with some Dirichlet boundaries

and test the value of the bottom-right corner cell.

>>> print numerix.allclose(phi(((L,), (0,))), valueBottomRight, atol = 1e-2)
1
>>> if __name__ == '__main__':
...     raw_input("Implicit steady-state diffusion. Press <return> to proceed...")