On Jul 18, 2008, at 6:59 AM, Etienne Rivard wrote:
Jonathan, although I'm starting a new thread here, I first want to
thank you for your once again very complete answer on the use of
operator variables to get access to additional functions. I'm not
100\% sure yet whether I'll need to use one of those. Nevertheless,
I found your answer extremely useful, it helped me a lot to
understand more of the mechanics of FiPy.
I'm glad it was helpful. Now we just need to figure out how to
properly document it and make it a bit more transparent for users to
actually do.
Now, I'd like to have your thoughts on another aspect of my problem
which is basically a Stefan problem or a free or moving boundary
problem, whatever you like to call it. I want to use a front-fixing
method to alleviate the moving boundary difficulty.
We wrote FiPy for solving phase field and level set problems, because
we didn't ever want to think about "the moving boundary difficulty"! 8^)
Here is an example, not the problem I want to solve, but the idea is
the same:
\[
\frac{\partial^2 u}{\partial \xi ^2} = s^2 \frac{\partial u}
{\partial t}
- s \xi \frac{ds}{dt}\frac{\partial u}{\partial \xi}
\]
I believe you'll want to divide through by $s^2$ and render this as
\[
\nabla\cdot\left(\frac{1}{s^2} u \right) = \frac{\partial u}{\partial t}
- \nabla\cdot\left(\xi \frac{1}{s} \frac{d s}{d t} u \right)
+ \frac{1}{s} \frac{d s}{d t} u
\]
or
DiffusionTerm(coeff=1/s**2) == TransientTerm() -
PowerLawConvectionTerm(coeff=xi * sdot / s) +
ImplicitSourceTerm(coeff=sdot / s)
\[
-\frac{1}{s} \left. \frac{\partial u}{\partial \xi}\right|_{\xi = 1}
= \frac{ds}{dt}
\]
You want to use this as an ODE to solve for s or as a boundary
condition on u? Either way, you need at least one more BC on u.
My first idea was to declare $s$ as a Variable() as in the diffusion
example on page 64 of the manual, ``manually'' linearize $\frac{ds}
{dt}$ with something like $\frac{s-s_{old}}{\Delta t}$ and solve the
resulting algebraic equation for $s$ at every sweep.
You'll need to explicitly track $s_{old}$, as only CellVariables have
that property. Alternatively, you could set up a one-element mesh and
make $s$ a CellVariable on that mesh and solve it as a PDE, but we're
pretty sure you'll have trouble communicating values back and forth
between the two different meshes. I think FiPy will raise an error if
you try to use $s$, defined on one mesh, as a coefficient for solving
a variable on a different mesh, but maybe you'll get lucky... no,
unfortunately you won't. I just tried:
>>> from fipy import *
>>> mesh = Grid1D(nx=10)
>>> xi = mesh.getCellCenters()
>>> u = CellVariable(mesh=mesh, hasOld=True)
>>> smesh = Grid1D(nx=1)
>>> s = CellVariable(mesh=smesh, value=1. hasOld=True)
>>> dt = 1.
>>> sdot = (s - s.getOld()) / dt
>>> eq = DiffusionTerm(coeff=1/s**2) == TransientTerm() -
PowerLawConvectionTerm(coeff=xi * sdot / s) +
ImplicitSourceTerm(coeff=sdot / s)
I was shocked to actually get this far, but when I tried:
>>> s.updateOld()
>>> s.setValue(2.)
>>> eq.solve(vars=u)
it throws "ValueError: shape mismatch: objects cannot be broadcast to
a single shape"
It turns out this is because both DiffusionTerm and ConvectionTerm
want their coefficients as FaceVariables, and will interpolate a
CellVariable to a FaceVariable if necessary, but a FaceVariable on
smesh has two faces, instead of the optimal 1 or acceptable 11. I can
probably come up with an ugly hack to avoid this, but will need to
think about it.
To get $\left. \frac{\partial u}{\partial \xi}\right|_{\xi = 1}$, I
think you may need to resort to the techniques discussed in this
recent thread:
http://thread.gmane.org/gmane.comp.python.fipy/973
but maybe Daniel or Recif has something better?
