Hello Jonathan,
Thank you very much for - once again - giving such a complete answer. I
think it's just fantastic that you guys care so much.
Unfortunately, for several reasons I have decided to go with a different
approach. Because I have a moving boundary problem, I have to modify the
equations to account for that. And after this coordinate transformation
and the one used to account for spherical symmetry, the equations change
considerably from their original state of convection-diffusion.
I will use the "method of lines". It's fairly easy to implement in 1-D
and allows the use of a generic ODE solver. For some reason it also
gives me a feeling that I'm grasping more of what's going on.
I'm not forgetting about FiPy though and it's quite possible that I need
it in the future. Thank you for your precious help.
Etienne Rivard
Jonathan Guyer wrote:
On Feb 26, 2008, at 7:16 AM, Etienne Rivard wrote:
After looking at the classes of various terms in the examples, I
couldn't find a way to represent $\frac{\partial}{\partial \xi}$. The
module examples.phase.impingement.mesh40x1.input has a similar term
but it uses the absolute, or magnitude of theta and this is not
exactly what I have. I have tried removing the getMag() part so that I
define theta.getGrad(), but that doesn't work, it's not the right class.
It's clearer (to me, at least) to represent the spatial derivatives
consistently, either as $\frac{\partial}{\partial \xi}$ or as $\nabla$,
but not mixed. My preference is for $\nabla$, because it highlights the
vector relationships that are tripping you up, I think. E.g.,
\[ \xi^2 \frac{\partial}{\partial t} (cy_F) - \xi^3 \frac{\dot{R}}{R}
\cdot \nabla (cy_F) + \frac{1}{R} \nabla \cdot
(\xi^2 c \vec{v} y_F) = \frac{1}{R^2} \nabla \cdot (\xi^2 c \bar{D}
\nabla y_F) \]
\[ \xi^2 \frac{\partial}{\partial t} (cy_F \theta) - \xi^3
\frac{\dot{R}}{R} \cdot \nabla (cy_F \theta) +
\frac{1}{R} \nabla \cdot (\xi^2 c \vec{v} y_F \theta) = \frac{1}{R^2}
\nabla \cdot (\xi^2 c \tilde{D} y_F \nabla \theta) + \frac{1}{R^2}
\nabla \cdot (\xi^2 c \tilde{D} \theta \nabla y_F) \]
\[ \xi^2 \frac{\partial}{\partial t} (cy_F \Psi) - \xi^3
\frac{\dot{R}}{R} \cdot \nabla (cy_F \Psi) +
\frac{1}{R} \nabla \cdot (\xi^2 c \vec{v} y_F \Psi) = \frac{1}{R^2}
\nabla
\cdot (\xi^2 c \hat{D} y_F \nabla \Psi) + \frac{1}{R^2} \nabla \cdot
(\xi^2 c \hat{D} \Psi \nabla y_F) \]
so that you can see that the second term, resulting from the material
derivative, is a velocity vector $\xi \frac{\dot{R}}{R}$ dotted with a
field gradient $\nabla (cy_F)$, etc. When you multiplied
theta.getGrad() times the scalar quantity xi**3 * Rdot / R, you got a
vector source, which is incompatible with your scalar governing
equation. As far as FiPy is concerned, even in 1D, a scalar field and a
vector field are different, even if they hold the same values. So, the
first step is to make a vector variable to hold your velocity, and then
dot that with the gradient. Assuming you're using FiPy 1.2, the simplest
(but not best) way is:
xi = mesh.getCellCenters()[...,0]
velocity = VectorCellVariable(mesh=mesh)
velocity[...,0] = xi * Rdot / R
and then represent the second term of the first equation as
xi**2 * velocity.dot(c * yF.getGrad())
I didn't set the "velocity" to xi**3 * Rdot / R because it isn't
dimensionally accurate, but you might find that it makes things a little
easier.
Although the dot form will work, we don't recommend it. It's a
first-order difference and, depending on magnitude, it can destabilize
your equations. In my own work, I've made the transformation $\vec{v}
\cdot \nabla \phi \rightarrow \nabla \cdot \vec{v} \phi + \phi \cdot
\nabla \cdot \vec{v}$, which allows you to use an implicit convection
term and an implicit source term, instead of an explicit source. To do
this, you'll need to evaluate "velocity" on the faces:
xiF = mesh.getFaceCenters()[...,0]
velocityF = VectorFaceVariable(mesh=mesh)
velocityF[...,0] = xiF * Rdot / R
and then represent the second term as
PowerLawConvectionTerm(coeff=c * xiF**2 * velocity) -
ImplicitSourceTerm(coeff=c * (xiF**2 * velocity).getDivergence())
Please let me know if I've misunderstood what you're trying to do, or if
this doesn't work for you.
If you want to transform your results in $\xi$ back into the $r$ frame,
we did some similar transformations between laboratory and reference
frames in a recent paper
@Article{Boettinger:2007p33,
author = {W. J Boettinger and J. E Guyer and C. E Campbell and G. B
McFadden},
journal = {P R Soc A},
pages = {3347--3373},
volume = {463},
year = {2007},
month = {Jan},
language = {English},
keywords = {Similarity Solution, Diffusion, Plane, Deformation,
Interdiffusion, Migration, Boltzmann-Matano, Phase Field, Dictra},
doi = {10.1098/rspa.2007.1904},
URL = {http://journals.royalsociety.org/content/n37t3346g7g17787/},
}
We included the heavily-commented FiPy script that we used for the phase
field portion in the electronic supplementary material.
Let me know if you can't obtain it and I'll send you a copy.
