Hi all,I'm new to FiPy and FVA in general, and am excited to use it to solve some thermal diffusion problems in my work. I figured before I tackled a full-scale problem, I'd start small. Unfortunately I was not very successful in this, and ran into an unexpected problem with the solution diverging.
The problem I set up is the diffusion equation for a simple 3D cylinder with square cross section and an aspect ratio of 2:1:1, with BCs of -1 and 1 on the two opposite long ends.
This works fine using FiPy's built-in Grid3D (e.g. attached script, fipytest5.py), and the solution happily approaches the expected steady state. However, I will need to solve problems with more complicated geometries, for which I expect to use gmsh. To test gmsh, I created the same shape in gmsh (fipytest6.geo), let gmsh handle the meshing, and let FiPy run on this mesh (fipytest6.py), expecting essentially the same solution.
However, it's clear that the gmsh solution starts to diverge after 60 or so steps and continuing much beyond that gives a completely nonsensical result. Reducing the step size does not help. Increasing the mesh density makes the problem worse (it diverges after fewer steps).
I've run the 2D version of this problem with both Grid2D and gmsh meshes (in analogy with circle.py), and they work fine, so I seem to be screwing up the transition to 3D.
Any insight would be greatly appreciated. Versions: In [283]: print sys.version 2.6.4 (r264:75706, Jun 4 2010, 18:20:31) [GCC 4.4.4 20100503 (Red Hat 4.4.4-2)] In [284]: print fipy.__version__ 2.1.2 [~]$ uname -r 2.6.34.9-69.fc13.x86_64 [~]$ gmsh --version 2.5.0 Thanks, Justin
#!/bin/env python
from fipy import *
dx = dy = dz = 1.
nz = ny = 20
nx = 40
lc = 1.
mesh5 = Grid3D(dx=dx, dy=dy, dz=dz, nx=nx, ny=ny, nz=nz)
phi5 = CellVariable(name='\phi_5', mesh=mesh5, value=0.)
eq5 = TransientTerm() == DiffusionTerm()
leftBC, rightBC = -1., 1.
# BCs5 = (FixedValue(faces=mesh5.getFacesLeft(), value=leftBC),
# FixedValue(faces=mesh5.getFacesRight(), value=rightBC))
x, y, z = mesh5.getFaceCenters()
leftFaces = mesh5.getExteriorFaces() & (x < min(x) + 0.05)
rightFaces = mesh5.getExteriorFaces() & (x > max(x) - 0.05)
BCs5 = (FixedValue(faces=leftFaces, value=leftBC),
FixedValue(faces=rightFaces, value=rightBC))
dt5 = 0.1
steps5 = 55
def takestep(phi, eq, BCs, dt):
eq.solve(var=phi, boundaryConditions=BCs, dt=dt)
def ts5(num=0):
"""Advance \phi_5 by num steps."""
for i in range(num):
takestep(phi5, eq5, BCs5, dt5)
ts5(steps5)
# for step in range(steps):
# takestep(phi5, eq, BCs, dt)
# # eq.solve(var=phi5, boundaryConditions=BCs, dt=timeStepDuration)
viewer = None
def view(phi):
"""Make a viewer to view the mesh phi."""
global viewer
viewer = Viewer(vars=phi, datamin=-1., datamax=1.)
fipytest6.geo
Description: application/vnd.dynageo
#!/bin/env python
from fipy import *
mesh6 = GmshImporter3D('fipytest6.msh')
phi6 = CellVariable(name='\phi_6', mesh=mesh6, value=0.)
eq6 = TransientTerm() == DiffusionTerm()
leftBC, rightBC = -1., 1.
# Doesn't seem to set one of the side's BCs correctly...
# BCs6 = (FixedValue(faces=mesh6.getFacesLeft(), value=leftBC),
# FixedValue(faces=mesh6.getFacesRight(), value=rightBC))
# But this works.
x, y, z = mesh6.getFaceCenters()
leftFaces = mesh6.getExteriorFaces() & (x < min(x) + 0.05)
rightFaces = mesh6.getExteriorFaces() & (x > max(x) - 0.05)
BCs6 = (FixedValue(faces=leftFaces, value=leftBC),
FixedValue(faces=rightFaces, value=rightBC))
dt6 = 0.01
steps6 = 550
def takestep(phi, eq, BCs, dt):
eq.solve(var=phi, boundaryConditions=BCs, dt=dt)
def ts6(num=0):
"""Advance \phi_6 by num steps."""
for i in range(num):
takestep(phi6, eq6, BCs6, dt6)
ts6(steps6)
viewer = None
def view(phi):
"""Make a viewer to view the mesh phi."""
global viewer
viewer = Viewer(vars=phi, datamin=-1., datamax=1.)
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
