Jonathan Guyer wrote:
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.
I'll trust you on that ;-)
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^)
Right. I remember once considering swithcing to another method like the
phase field method to solve my problem. But I have to admit I was
somewhat demotivated by the apparent complexity of the method. For the
moment, I don't think I'll get deeper into that.
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.
Thanks for the offer, I appreciate it but. I have time to work on this
and this actually a part of my PhD work so I need to show a deep
understanding of everything I use, so I would rather shoot for a more
elegant solution.
Also, it makes me feel really good to know that this problem is not as
simple as it looks and actually requires some thought.
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?
I'll have a close look at that. Thanks for the valuable input.
Etienne
