On Apr 10, 2007, at 1:33 PM, Elizabeth Tseng wrote:
I'm new to both fipy and PDEs.
Welcome to both.
I'm trying to construct a transient term that looks like this:
\frac{\partial PHI}{\partial t} = (D / (r^2)) \frac{r^2 \frac
{\partial PHI}{\partial r}}{\partial r}
One minor correction:
\frac{\partial PHI}{\partial t} = (D / (r^2)) \frac{\partial r^2 \frac
{\partial PHI}{\partial r}}{\partial r}
^^^^^^^^
or
\frac{\partial PHI}{\partial t} = (D / (r^2)) \frac{\partial}
{\partial r} (r^2 \frac{\partial PHI}{\partial r})
I know the left hand term should be a Transient Term. But what
about the right hand term? The coefficient, instead of being a
fixed constant D, is D / (r^2)...
The right hand term is a Laplacian of D PHI in spherical coordinates,
so this is still a DiffusionTerm, but capturing the radial dependence
is a little tricky. One way to deal with this in FiPy is with a
conical or wedge-shaped mesh that gets wider as you move away from
the origin. This is something that we've talked about before and we
should put one together just to see how it performs.
We realized, however, that there's an easier way to deal with this
when all you have is radial dependence. You can multiply both sides
of the equation by r^2, and noting that there is no time-dependence
to r and no r-dependence to D, we can write:
\frac{\partial r^2 PHI}{\partial t} = \frac{\partial}{\partial r} (D
r^2 \frac{\partial PHI}{\partial r})
so we have a TransientTerm with coefficient r^2 on the left and a
DiffusionTerm with a coefficient of D r^2 on the right. Here's a way
to do that in FiPy:
>>> from fipy import *
>>> mesh = Grid1D(nx=100, dx=1)
>>> phi = CellVariable(mesh=mesh)
>>> r = mesh.getCellCenters()[...,0]
>>> r2 = CellVariable(mesh=mesh, value=r**2)
>>> phi.setValue(1, where=r < 50)
>>> D = 1
>>> eq = TransientTerm(coeff=r2) == DiffusionTerm(coeff=D * r2)
>>> viewer = viewers.make(vars=phi, limits={'datamin':0,
'datamax':1})
>>> viewer.plot()
>>> for step in range(1000):
... eq.solve(var=phi)
... viewer.plot()
If you need to worry about the \phi and \theta dependences, too, then
this rearrangement won't work, but for purely radial problems, it
should be fine.
