Merge branch 'ticket497-examples_phase_diffusion_py_has_problems' int…

…o develop

Fixes ticket:497

examples/phase/binary.py

@@ -107,7 +107,7 @@
 variety of choices for the interpolation function :math:`p(\phi)` and the
 barrier function :math:`g(\phi)`, 
    
-such as those shown in mod:`examples.phase.simple`
+such as those shown in :mod:`examples.phase.simple`
 
 >>> def p(phi):
 ...     return phi**3 * (6 * phi**2 - 15 * phi + 10)
@@ -374,7 +374,7 @@ def deltaChemPot(phase, C, T):
 
 >>> Dl = Variable(value=1e-5) # cm**2 / s
 >>> Ds = Variable(value=1e-9) # cm**2 / s
->>> D = (Dl - Ds) * phase.arithmeticFaceValue + Dl
+>>> D = (Ds - Dl) * phase.arithmeticFaceValue + Dl
 
 >>> phaseTransformationVelocity = \
 ...  ((enthalpyB - enthalpyA) * p(phase).faceGrad
@@ -508,7 +508,7 @@ def deltaChemPot(phase, C, T):
 earlier examples that the diffusion problem is unconditionally stable, we
 need take only one very large timestep to reach equilibrium
 
->>> dt = 1.e2
+>>> dt = 1.e5
 
 Because the phase field equation is coupled to the composition through
 ``enthalpy`` and ``W`` and the diffusion equation is coupled to the phase
@@ -575,7 +575,35 @@ def deltaChemPot(phase, C, T):
 (solidify), we will need to take much smaller timesteps (the time scales of
 diffusion and of phase transformation compete with each other).
 
->>> dt = 1.e-6
+The `CFL limit`_ requires that no interface should advect more than one grid
+spacing in a timestep. We can get a rough idea for the maximum timestep we can
+take by looking at the velocity of convection induced by phase transformation in
+Eq. :eq:`eq:phase:binary:diffusion:canonical`. If we assume that the phase changes from 1 to 0 in a single grid spacing,
+that the diffusivity is `Dl` at the interface, and that the term due to the difference in
+barrier heights is negligible:
+    
+.. math::
+
+   \vec{u}_\phi &= \frac{D_\phi}{C} \nabla \phi
+   \\
+   &\approx
+   \frac{Dl \frac{1}{2} V_m}{R T}
+   \left[
+       \frac{L_B\left(T - T_M^B\right)}{T_M^B} 
+       - \frac{L_A\left(T - T_M^A\right)}{T_M^A}
+   \right] \frac{1}{\Delta x}
+   \\
+   &\approx
+   \frac{Dl \frac{1}{2} V_m}{R T}
+   \left(L_B + L_A\right) \frac{T_M^A - T_M^B}{T_M^A + T_M^B} 
+   \frac{1}{\Delta x}
+   \\
+   &\approx \unit{0.28}{\centi\meter\per\second}
+
+To get a :math:`\text{CFL} = \vec{u}_\phi \Delta t / \Delta x < 1`, we need a
+time step of about :math:`\unit{10^{-5}}{\second}`.
+
+>>> dt = 1.e-5
 
 >>> if __name__ == '__main__':
 ...     timesteps = 100
@@ -618,8 +646,7 @@ def deltaChemPot(phase, C, T):
    examine different temperatures in this example, so we declare :math:`T` 
    as a :class:`~fipy.variables.variable.Variable`
 
-.. .. bibmissing:: /documentation/refs.bib
-    :sort:
+.. _CFL limit: http://en.wikipedia.org/wiki/Courant-Friedrichs-Lewy_condition
 """
 
 __docformat__ = 'restructuredtext'

examples/phase/binaryCoupled.py

@@ -108,7 +108,7 @@
 variety of choices for the interpolation function :math:`p(\phi)` and the
 barrier function :math:`g(\phi)`, 
    
-such as those shown in mod:`examples.phase.simple`
+such as those shown in :mod:`examples.phase.simple`
 
 >>> def p(phi):
 ...     return phi**3 * (6 * phi**2 - 15 * phi + 10)
@@ -385,7 +385,7 @@ def deltaChemPot(phase, C, T):
 
 >>> Dl = Variable(value=1e-5) # cm**2 / s
 >>> Ds = Variable(value=1e-9) # cm**2 / s
->>> Dc = (Dl - Ds) * phase.arithmeticFaceValue + Dl
+>>> Dc = (Ds - Dl) * phase.arithmeticFaceValue + Dl
 
 >>> Dphi = ((Dc * C.harmonicFaceValue * (1 - C.harmonicFaceValue) * Vm / (R * T))
 ...         * ((enthalpyB - enthalpyA) * pPrime(phase.arithmeticFaceValue)
@@ -517,7 +517,7 @@ def deltaChemPot(phase, C, T):
 earlier examples that the diffusion problem is unconditionally stable, we
 need take only one very large timestep to reach equilibrium
 
->>> dt = 1.e2
+>>> dt = 1.e5
 
 Because the phase field equation is coupled to the composition through
 ``enthalpy`` and ``W`` and the diffusion equation is coupled to the phase
@@ -587,7 +587,37 @@ def deltaChemPot(phase, C, T):
 (solidify), we will need to take much smaller timesteps (the time scales of
 diffusion and of phase transformation compete with each other).
 
->>> dt = 1.e-6
+The `CFL limit`_ requires that no interface should advect more than one grid
+spacing in a timestep. We can get a rough idea for the maximum timestep we can
+take by looking at the velocity of convection induced by phase transformation in
+Eq. :eq:`eq:phase:binary:diffusion:canonical` (even though there is no explicit
+convection in the coupled form used for this example, the principle remains the
+same). If we assume that the phase changes from 1 to 0 in a single grid spacing,
+that the diffusivity is `Dl` at the interface, and that the term due to the difference in
+barrier heights is negligible:
+    
+.. math::
+
+   \vec{u}_\phi &= \frac{D_\phi}{C} \nabla \phi
+   \\
+   &\approx
+   \frac{Dl \frac{1}{2} V_m}{R T}
+   \left[
+       \frac{L_B\left(T - T_M^B\right)}{T_M^B} 
+       - \frac{L_A\left(T - T_M^A\right)}{T_M^A}
+   \right] \frac{1}{\Delta x}
+   \\
+   &\approx
+   \frac{Dl \frac{1}{2} V_m}{R T}
+   \left(L_B + L_A\right) \frac{T_M^A - T_M^B}{T_M^A + T_M^B} 
+   \frac{1}{\Delta x}
+   \\
+   &\approx \unit{0.28}{\centi\meter\per\second}
+
+To get a :math:`\text{CFL} = \vec{u}_\phi \Delta t / \Delta x < 1`, we need a
+time step of about :math:`\unit{10^{-5}}{\second}`.
+
+>>> dt = 1.e-5
 
 >>> if __name__ == '__main__': 
 ...     timesteps = 100
@@ -627,7 +657,7 @@ def deltaChemPot(phase, C, T):
    examine different temperatures in this example, so we declare :math:`T` 
    as a :class:`~fipy.variables.variable.Variable`
 
-
+.. _CFL limit: http://en.wikipedia.org/wiki/Courant-Friedrichs-Lewy_condition
 """
 
 __docformat__ = 'restructuredtext'