examples.diffusion.circle

Solve the diffusion equation in a circular domain meshed with triangles.

This example demonstrates how to solve a simple diffusion problem on a non-standard mesh with varying boundary conditions. The Gmsh package is used to create the mesh. Firstly, define some parameters for the creation of the mesh,

>>> cellSize = 0.05
>>> radius = 1.

The cellSize is the preferred edge length of each mesh element and the radius is the radius of the circular mesh domain. In the following code section a file is created with the geometry that describes the mesh. For details of how to write such geometry files for Gmsh, see the gmsh manual.

The mesh created by Gmsh is then imported into FiPy using the Gmsh2D object.

>>> from fipy import *
>>> mesh = Gmsh2D('''
...               cellSize = %(cellSize)g;
...               radius = %(radius)g;
...               Point(1) = {0, 0, 0, cellSize};
...               Point(2) = {-radius, 0, 0, cellSize};
...               Point(3) = {0, radius, 0, cellSize};
...               Point(4) = {radius, 0, 0, cellSize};
...               Point(5) = {0, -radius, 0, cellSize};
...               Circle(6) = {2, 1, 3};
...               Circle(7) = {3, 1, 4};
...               Circle(8) = {4, 1, 5};
...               Circle(9) = {5, 1, 2};
...               Line Loop(10) = {6, 7, 8, 9};
...               Plane Surface(11) = {10};
...               ''' % locals()) 

Using this mesh, we can construct a solution variable

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

We can now create a Viewer to see the mesh

>>> viewer = None
>>> if __name__ == '__main__':
...     try:
...         viewer = Viewer(vars=phi, datamin=-1, datamax=1.)
...         viewer.plotMesh()
...         raw_input("Irregular circular mesh. Press <return> to proceed...") 
...     except:
...         print "Unable to create a viewer for an irregular mesh (try Gist2DViewer, Matplotlib2DViewer, or MayaviViewer)"

We set up a transient diffusion equation

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

The following line extracts the coordinate values on the exterior faces. These are used as the boundary condition fixed values.

>>> X, Y = mesh.faceCenters 

>>> phi.constrain(X, mesh.exteriorFaces) 

We first step through the transient problem

>>> timeStepDuration = 10 * 0.9 * cellSize**2 / (2 * D)
>>> steps = 10
>>> for step in range(steps):
...     eq.solve(var=phi,
...              dt=timeStepDuration) 
...     if viewer is not None:
...         viewer.plot() 

---

If we wanted to plot or analyze the results of this calculation with another application, we could export tab-separated-values with

TSVViewer(vars=(phi, phi.grad)).plot(filename="myTSV.tsv")

x	y	solution variable	solution variable_grad_x	solution variable_grad_y
0.975559734792414	0.0755414402612554	0.964844363287199	-0.229687917881182	0.00757854476106331
0.0442864953037566	0.79191893162384	0.0375859836421991	-0.773936613923853	-0.205560697612547
0.0246775505084069	0.771959648896982	0.020853932412869	-0.723540342405813	-0.182589694334729
0.223345558247991	-0.807931073108895	0.203035857140125	-0.777466238738658	0.0401235242511506
-0.00726763301939488	-0.775978916110686	-0.00412895434496877	-0.650055516507232	-0.183112882869288
-0.0220279064527904	-0.187563765977912	-0.012771874945585	-0.35707168379437	-0.056072788439713
0.111223320911545	-0.679586798311355	0.0911595298310758	-0.613455176718145	0.0256182541329463
-0.78996770899909	-0.0173672729866294	-0.693887874335319	-1.00671109050419	-0.127611490372511
-0.703545986179876	-0.435813500559859	-0.635004192597412	-0.896203033957194	-0.00855563518923689
0.888641841567831	-0.408558914368324	0.877939107374768	-0.32195762184087	-0.22696791637322
0.38212257821916	-0.51732949653553	0.292889724306196	-0.854466141879776	0.199715815696975
-0.359068256998365	0.757882581524374	-0.323541041763627	-0.870534227755687	0.0792631912863636
-0.459673905457569	-0.701526587772079	-0.417577664032421	-0.725460726303266	-0.119132299176163
-0.338256179134518	-0.523565732643067	-0.254030052182524	-0.923505840608445	-0.192224240688976
0.87498754712638	0.174119064688993	0.836057900916614	-1.11590500805745	-0.211010116496191
-0.484106960369249	0.0705987421869745	-0.319827850867342	-0.867894407968447	0.051246727010685
-0.0221203060940465	-0.216026820080053	-0.0152729438559779	-0.341246696530392	-0.0538476142281317

The values are listed at the cell centers. Particularly for irregular meshes, no specific ordering should be relied upon. Vector quantities are listed in multiple columns, one for each mesh dimension.

---

This problem again has an analytical solution that depends on the error function, but it’s a bit more complicated due to the varying boundary conditions and the different horizontal diffusion length at different vertical positions

>>> x, y = mesh.cellCenters 
>>> t = timeStepDuration * steps

>>> phiAnalytical = CellVariable(name="analytical value",
...                              mesh=mesh) 

>>> x0 = radius * numerix.cos(numerix.arcsin(y)) 
>>> try:
...     from scipy.special import erf 
...     ## This function can sometimes throw nans on OS X
...     ## see http://projects.scipy.org/scipy/scipy/ticket/325
...     phiAnalytical.setValue(x0 * (erf((x0+x) / (2 * numerix.sqrt(D * t))) 
...                                  - erf((x0-x) / (2 * numerix.sqrt(D * t))))) 
... except ImportError:
...     print "The SciPy library is not available to test the solution to \
... the transient diffusion equation"

>>> print phi.allclose(phiAnalytical, atol = 7e-2) 
1

>>> if __name__ == '__main__':
...     raw_input("Transient diffusion. Press <return> to proceed...")

---

As in the earlier examples, we can also directly solve the steady-state diffusion problem.

>>> DiffusionTerm(coeff=D).solve(var=phi) 

The values at the elements should be equal to their x coordinate

>>> print phi.allclose(x, atol = 0.03) 
1

Display the results if run as a script.

>>> if viewer is not None:
...     viewer.plot()
...     raw_input("Steady-state diffusion. Press <return> to proceed...")