Hello again,
I am trying to solve the diffusion example in a circle with an annular region
instead. I thought this would be a straightforward extension and I’m sure it
is, I am just missing something. The reason I think something is wrong is
because the solutions are looking like this, after importing the data into
Mathematica:
The mesh I created seems to be ok however:
My code is below:
- I suspect that there is something wrong with the Boundary Conditions. Or
perhaps the “Physical” labels of the gmsh code. I saw something about this in
the Documentation, but I was unsure what was meant.
cellSize=0.05
radius=1.
from fipy import *
#Creates the annular mesh
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(5) = {-radius, 0, 0, cellSize};
Point(6) = {0, -radius, 0, cellSize};
Circle(1) = {2, 1, 3};
Circle(2) = {3, 1, 5};
Circle(3) = {5, 1, 6};
Circle(4) = {6, 1, 2};
Line Loop(1) = {1, 2, 3, 4};
Point(200) = {radius/2, 0, 0, cellSize};
Point(300) = {0, radius/2, 0, cellSize};
Point(500) = {-radius/2, 0, 0, cellSize};
Point(600) = {0, -radius/2, 0, cellSize};
Circle(100) = {200, 1, 300};
Circle(200) = {300, 1, 500};
Circle(300) = {500, 1, 600};
Circle(400) = {600, 1, 200};
Line Loop(100) = {100, 200, 300, 400};
Plane Surface(1)={1,100};
Physical Line(1001) = {1,2,3,4};
Physical Line(1002) = {100,200,300,400};
Physical Surface(1005) = {1};
''' % locals())
phi=CellVariable(name="solution",mesh=mesh,value=0.)
D=1
eq=TransientTerm()==DiffusionTerm(coeff=D)
X,Y=mesh.faceCenters
#Boundary Conditions
phi.constrain(X,mesh.exteriorFaces)
phi.constrain(X,mesh.interiorFaces)
timeStepDuration=10*0.9*cellSize**2/(2*D)
steps=10
for step in range(steps):
eq.solve(var=phi,dt=timeStepDuration)
TSVViewer(vars=(phi)).plot(filename="diffusion.annulus."+str(step)+".tsv”)
Any and all help appreciated!
Thanks,
Kyle_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
