Objectives:
First thing we'll do is use the qt backend for interacting with matplotlib. There is documentation in matplotlib about backends. One or more of the available backends may be installed with your python distribution.
In [1]:
%matplotlib osx
from fipy import *
In [ ]:
%matplotlib
from fipy import *
In [2]:
L = 1.
nx = 200
dx = L / nx
timeStepDuration = 0.001
steps = 100
Here we set the diffusion coefficient.
In [3]:
D11 = 0.5
Note: The 'mesh' command creates the mesh (gridpoints) on which we will solve the equation. This is specific to FiPy.
At this point, if you are in the IPython notebook I would suggest you try the following in the cell below:
This is a powerful way to explore the available functions.
In [4]:
mesh = Grid1D(dx = dx, nx = nx)
We assign to 'c' objects that are 'CellVariables'. This is a special type of variable used by FiPy to hold the values for the concentration in our solution. We also create a viewer here so that we can inspect the values of 'c'.
In [5]:
c = CellVariable(name = "c", mesh = mesh)
viewer = MatplotlibViewer(vars=(c,),datamin=-0.1, datamax=1.1, legend=None)
Note: This command sets 'x' to contain a list of numbers that define the x position of the grid-points.
In [6]:
x = mesh.cellCenters
x
Out[6]:
In [ ]:
x
In [7]:
c.setValue(0.0)
viewer.plot()
In [35]:
c.setValue(0.2, where=x < L/2.)
c.setValue(0.8, where=x > L/2.)
viewer.plot()
Note: Boundary conditions come in two types. Fixed flux and fixed value. The syntax is:
Fixed value boundaries, can set VALUE or FLUX as a float
In [15]:
boundaryConditions=(FixedValue(mesh.facesLeft,0.2),
FixedValue(mesh.facesRight,0.8))
This line defines the diffusion equation:
In [16]:
eqn1 = TransientTerm() == ImplicitDiffusionTerm(D11)
In [36]:
for step in range(1000):
eqn1.solve(c, boundaryConditions = boundaryConditions, dt = timeStepDuration)
viewer.plot()
These lines print out a text file with the final data
In [37]:
from fipy.viewers.tsvViewer import TSVViewer
TSVViewer(vars=(c)).plot(filename="output.txt")
In [38]:
!head output.txt
In [ ]:
%matplotlib osx
from fipy import *
nx is the number of grid points we wish to have in our solution. If we have fewer we sacrifice accuracy, if we have more, the computational time increases and we are subject to roundoff error. You should always check that your solution does not depend on the number of grid points and the grid spacing!dx is the spacing between grid points. Similar comments as above.Accuracy and stability of the solution depend on choices for the timeStepDuration and the grid point spacing and the diffusion coefficient. We will not discuss stability any further, just know that it is something that needs to be considered.
In [39]:
L = 1.
nx = 50
dx = L / nx
timeStepDuration = 0.001
steps = 100
Note: In the first equation, the diffusion coefficient is constant, concentration independent.
In [51]:
D1 = 3.0
Note: You have seen all of the following code before. This time we are solving two simultaneous equations, eqn1 and eqn2.
Note: The 'mesh' command creates the mesh (gridpoints) on which we will solve the equation. This is specific to FiPy.
In [41]:
mesh = Grid1D(dx = dx, nx = nx)
We assign to 'c' objects that are 'CellVariables'. This is a special type of variable used by FiPy to hold the values for the concentration in our solution. c1 for eqn1 and c2 for eqn2
In [42]:
c1 = CellVariable(
name = "c1",
mesh = mesh,
hasOld = True)
c2 = CellVariable(
name = "c2",
mesh = mesh,
hasOld = True)
Note: This command sets 'x' to contain a list of numbers that define the x position of the grid-points.
In [43]:
x = mesh.cellCenters
These lines set diffusant in the initial condition. They are set to the same values for easy comparison. Feel free to change these values.
In [55]:
c1.setValue(0.8)
c1.setValue(0.2, where=x > L/3.)
c2.setValue(0.8)
c2.setValue(0.2, where=x > L/3.)
In [56]:
viewer.plot()
In [53]:
viewer = Matplotlib1DViewer(vars = (c1,c2), limits = {'datamin': 0., 'datamax': 1.})
viewer.plot()
Boundary conditions can be either fixed flux or fixed value. Here, fixed value is used for simple comparison between the diffusion profiles.
In [46]:
boundaryConditions=(FixedValue(mesh.facesLeft,0.8), FixedValue(mesh.facesRight,0.2))
In [47]:
boundaryConditions=(FixedFlux(mesh.facesLeft,0.0), FixedFlux(mesh.facesRight,0.0))
Note: In the second equation, the diffusion coefficient is non-constant and is a function of concentration in the system. So we use D22_0 and D22_1 as the end points of our function. The function, given by "D" is simply a linear interpolation of the two D values. You are, of course, free to try other functions.
In [48]:
D22_0 = 3.0
D22_1 = 0.5
D2 = (D22_1 - D22_0)*c2 + D22_0
Note: These are the two diffusion equations. The first equation is as in previous code, using a concentration independent D1. The second equation uses a non-constant D described above.
In [49]:
eqn1 = TransientTerm() == ImplicitDiffusionTerm(D1)
eqn2 = TransientTerm() == ImplicitDiffusionTerm(D2)
The following is an if loop that waits for user input before executing. Iterates for number of steps stated earlier for each equation.
In [57]:
for step in range(10):
c1.updateOld()
c2.updateOld()
eqn1.solve(c1, boundaryConditions = boundaryConditions, dt = timeStepDuration)
eqn2.solve(c2, boundaryConditions = boundaryConditions, dt = timeStepDuration)
viewer.plot()
In [ ]: