Re: Front-fixing method for moving boundary problem
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
RSS Feed