From f0ce533ba4cdccfc68c642449f393977733feb0f Mon Sep 17 00:00:00 2001 From: Jonathan Guyer Date: Thu, 12 Sep 2013 10:01:17 -0400 Subject: [PATCH 1/3] Fix phase-weighted diffusivity and increase time steps Diffusivity did not approach proper bulk limits. Correct weighting allows larger timesteps. Addresses ticket:497 --- examples/phase/binary.py | 6 +++--- examples/phase/binaryCoupled.py | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/phase/binary.py b/examples/phase/binary.py index e5cae032fc..8d78afc0a3 100755 --- a/examples/phase/binary.py +++ b/examples/phase/binary.py @@ -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,7 @@ 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 +>>> dt = 1.e-5 >>> if __name__ == '__main__': ... timesteps = 100 diff --git a/examples/phase/binaryCoupled.py b/examples/phase/binaryCoupled.py index 6c501e4731..8f1245c6ef 100755 --- a/examples/phase/binaryCoupled.py +++ b/examples/phase/binaryCoupled.py @@ -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,7 @@ 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 +>>> dt = 1.e-5 >>> if __name__ == '__main__': ... timesteps = 100 From 5406e13684e7d97a1148e6d652c598478ab809b5 Mon Sep 17 00:00:00 2001 From: Jonathan Guyer Date: Thu, 12 Sep 2013 10:08:21 -0400 Subject: [PATCH 2/3] Fix typo --- examples/phase/binary.py | 2 +- examples/phase/binaryCoupled.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/phase/binary.py b/examples/phase/binary.py index 8d78afc0a3..a7d8a0e613 100755 --- a/examples/phase/binary.py +++ b/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) diff --git a/examples/phase/binaryCoupled.py b/examples/phase/binaryCoupled.py index 8f1245c6ef..48d40018ba 100755 --- a/examples/phase/binaryCoupled.py +++ b/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) From fc62e8f3896ff20d11c9cbdd2b70a796d5bb4561 Mon Sep 17 00:00:00 2001 From: Jonathan Guyer Date: Fri, 13 Sep 2013 09:20:50 -0400 Subject: [PATCH 3/3] Document calculation of CFL Addresses ticket:497 --- examples/phase/binary.py | 31 +++++++++++++++++++++++++++++-- examples/phase/binaryCoupled.py | 32 +++++++++++++++++++++++++++++++- 2 files changed, 60 insertions(+), 3 deletions(-) diff --git a/examples/phase/binary.py b/examples/phase/binary.py index a7d8a0e613..1ba636c2c7 100755 --- a/examples/phase/binary.py +++ b/examples/phase/binary.py @@ -575,6 +575,34 @@ 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). +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__': @@ -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' diff --git a/examples/phase/binaryCoupled.py b/examples/phase/binaryCoupled.py index 48d40018ba..1388a680a6 100755 --- a/examples/phase/binaryCoupled.py +++ b/examples/phase/binaryCoupled.py @@ -587,6 +587,36 @@ 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). +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__': @@ -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'