Re: Front-fixing method for moving boundary problem

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?