Commit: 3c24066b5cd0bb10668cef2bf5c45939cd99b34e - myprojects-hg-reborn (hg)

Commit MetaInfo

• Revision: 3c24066b5cd0bb10668cef2bf5c45939cd99b34e (tree)
• Time: 2011-03-29 01:25:56
• Author: lorenzo
• Commiter: lorenzo

Log Message

I added a script solving the diffusion equation on a circle.

Change Summary

• add: Python-codes/circle.py

Incremental Difference

diff -r e61011792ce9 -r 3c24066b5cd0 Python-codes/circle.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Python-codes/circle.py	Mon Mar 28 16:25:56 2011 +0000
@@ -0,0 +1,78 @@
+#!/usr/bin/env python
+
+
+from fipy import *
+
+
+cellSize = 0.05
+radius = 1.
+
+transient=0
+
+mesh = GmshImporter2D('''
+                       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())
+
+
+
+phi = CellVariable(name = "solution variable",mesh = mesh,value = 0.)
+
+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)"
+
+D = 1.
+
+X, Y = mesh.getFaceCenters()
+
+BCs = (FixedValue(faces=mesh.getExteriorFaces(), value=X),)
+
+
+if (transient==1):
+
+    eq = TransientTerm() == DiffusionTerm(coeff=D)
+
+
+    timeStepDuration = 10 * 0.9 * cellSize**2 / (2 * D)
+    steps = 10
+    for step in range(steps):
+        eq.solve(var=phi,boundaryConditions=BCs,dt=timeStepDuration)
+        if viewer is not None:
+            viewer.plot()
+
+
+    TSVViewer(vars=(phi, phi.getGrad())).plot(filename="myTSV.tsv")
+
+# SS solution
+
+
+DiffusionTerm(coeff=D).solve(var=phi,boundaryConditions=BCs)
+
+
+print "phi.allclose(X, atol = 0.02), ", phi.allclose(X, atol = 0.02)
+
+if viewer is not None:
+    viewer.plot()
+    raw_input("Steady-state diffusion. Press <return> to proceed...")
+
+
+
+
+print "So far so good"