-
Notifications
You must be signed in to change notification settings - Fork 146
/
circle.py
executable file
·182 lines (135 loc) · 5.85 KB
/
circle.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
r"""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 :term:`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 :term:`Gmsh`, see the `gmsh manual`_.
.. _gmsh manual: http://www.geuz.org/gmsh/doc/texinfo/gmsh.html
The mesh created by :term:`Gmsh` is then imported into :term:`FiPy` using the
:class:`~fipy.meshes.gmshMesh.Gmsh2D` object.
>>> from fipy import CellVariable, Gmsh2D, TransientTerm, DiffusionTerm, Viewer
>>> from fipy.tools import numerix
>>> 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()) # doctest: +GMSH
Using this mesh, we can construct a solution variable
.. index::
object: fipy.variables.cellVariable.CellVariable
>>> phi = CellVariable(name = "solution variable",
... mesh = mesh,
... value = 0.) # doctest: +GMSH
We can now create a :class:`Viewer <~fipy.viewers.viewer.AbstractViewer>` to see the mesh
>>> viewer = None
>>> from fipy import input
>>> if __name__ == '__main__':
... try:
... viewer = Viewer(vars=phi, datamin=-1, datamax=1.)
... viewer.plotMesh()
... input("Irregular circular mesh. Press <return> to proceed...") # doctest: +GMSH
... except:
... print("Unable to create a viewer for an irregular mesh (try Matplotlib2DViewer or MayaviViewer)")
.. image:: circleMesh.*
:width: 90%
:align: center
:alt: circular mesh generated by Gmsh
We set up a transient diffusion equation
.. index::
object: fipy.terms.transientTerm.TransientTerm
object: fipy.terms.implicitDiffusionTerm.DiffusionTerm
>>> D = 1.
>>> eq = TransientTerm() == DiffusionTerm(coeff=D)
The following line extracts the :math:`x` coordinate values on the exterior
faces. These are used as the boundary condition fixed values.
>>> X, Y = mesh.faceCenters # doctest: +GMSH
>>> phi.constrain(X, mesh.exteriorFaces) # doctest: +GMSH
We first step through the transient problem
>>> timeStepDuration = 10 * 0.9 * cellSize**2 / (2 * D)
>>> steps = 10
>>> from builtins import range
>>> for step in range(steps):
... eq.solve(var=phi,
... dt=timeStepDuration) # doctest: +GMSH
... if viewer is not None:
... viewer.plot() # doctest: +GMSH
.. image:: circleTransient.*
:width: 90%
:align: center
:alt: evolution of diffusion problem on a circular mesh
----
If we wanted to plot or analyze the results of this calculation with
another application, we could export tab-separated-values with
.. index::
object: fipy.viewers.tsvViewer.TSVViewer
::
TSVViewer(vars=(phi, phi.grad)).plot(filename="myTSV.tsv")
.. literalinclude:: myTSV.tsv
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 # doctest: +GMSH
>>> t = timeStepDuration * steps
>>> phiAnalytical = CellVariable(name="analytical value",
... mesh=mesh) # doctest: +GMSH
.. index::
module: scipy
single: sqrt; arcsin; cos
>>> x0 = radius * numerix.cos(numerix.arcsin(y)) # doctest: +GMSH
>>> try:
... from scipy.special import erf # doctest: +SCIPY
... ## 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))))) # doctest: +GMSH, +SCIPY
... 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)) # doctest: +GMSH, +SCIPY
1
>>> from fipy import input
>>> if __name__ == '__main__':
... 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) # doctest: +GMSH
The values at the elements should be equal to their `x` coordinate
>>> print(phi.allclose(x, atol = 0.03)) # doctest: +GMSH
1
Display the results if run as a script.
>>> from fipy import input
>>> if viewer is not None:
... viewer.plot()
... input("Steady-state diffusion. Press <return> to proceed...")
.. image:: circleSteadyState.*
:width: 90%
:align: center
:alt: steady-state solution to diffusion on a circular mesh
"""
from __future__ import unicode_literals
__docformat__ = 'restructuredtext'
if __name__ == '__main__':
import fipy.tests.doctestPlus
exec(fipy.tests.doctestPlus._getScript())