-
Notifications
You must be signed in to change notification settings - Fork 146
/
mesh20x20Coupled.py
executable file
·109 lines (77 loc) · 2.93 KB
/
mesh20x20Coupled.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
r"""Solve a coupled set of diffusion equations in two dimensions.
This example solves a diffusion problem and demonstrates the use of
applying boundary condition patches.
.. index:: Grid2D
>>> from fipy import CellVariable, Grid2D, Viewer, TransientTerm, DiffusionTerm
>>> from fipy.tools import numerix
>>> nx = 20
>>> ny = nx
>>> dx = 1.
>>> dy = dx
>>> L = dx * nx
>>> mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
We create a :class:`~fipy.variables.cellVariable.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(var=phi) == DiffusionTerm(coeff=D, var=phi)
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
.. index::
module: fipy.viewers
>>> 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
>>> from builtins import range
>>> for step in range(steps):
... eq.solve(dt=timeStepDuration)
... if __name__ == '__main__':
... viewer.plot()
.. image:: mesh20x20transient.*
:width: 90%
:align: center
We can test the value of the bottom-right corner cell.
>>> print(numerix.allclose(phi(((L,), (0,))), valueBottomRight, atol = 1e-2))
1
>>> from fipy import input
>>> if __name__ == '__main__':
... input("Implicit transient diffusion. Press <return> to proceed...")
----
We can also solve the steady-state problem directly
>>> DiffusionTerm(var=phi).solve()
>>> if __name__ == '__main__':
... viewer.plot()
.. image:: mesh20x20steadyState.*
:width: 90%
:align: center
and test the value of the bottom-right corner cell.
>>> print(numerix.allclose(phi(((L,), (0,))), valueBottomRight, atol = 1e-2))
1
>>> from fipy import input
>>> if __name__ == '__main__':
... input("Implicit steady-state diffusion. Press <return> to proceed...")
"""
from __future__ import unicode_literals
__docformat__ = 'restructuredtext'
##from fipy.tools.profiler.profiler import Profiler
##from fipy.tools.profiler.profiler import calibrate_profiler
if __name__ == '__main__':
import fipy.tests.doctestPlus
exec(fipy.tests.doctestPlus._getScript())