On Nov 9, 2012, at 8:42 PM, Adam Stone wrote: > I was reviewing the examples.phase.binaryCoupled example at > > http://www.ctcms.nist.gov/fipy/examples/phase/generated/examples.phase.binaryCoupled.html > > and while trying to understand the implicit treatment of dC/dt, I noticed > some possible problems. The first is that some of the methods in the > definitons Dphi and diffusionEq no longer seem to work (specifically > phase.arithmeticFaceValue, C.harmonicFaceValue, and passing a var argument to > the transient and diffusion terms). I noticed that the simpleBinary1D.py file > in my examples folder does work and has different versions of these two > definitions, but you may want to update the website to prevent confusion.
Actually, the reverse is true. That example is for FiPy 3.0 and your symptoms indicate you're still running FiPy 2.x. You need to upgrade. > My main question involves the definitions of diffusion coefficients which is > the same in both versions: > > Dl = Variable(value=1e-5) # cm**2 / s > Ds = Variable(value=1e-9) # cm**2 / s > D = (Dl - Ds) * phase.getArithmeticFaceValue() + Dl > > Perhaps I'm misunderstanding the basis for this expression, but to me it > seems like there is an error and that the correct expression should be > > D = (Ds - Dl) * phase.getArithmeticFaceValue() + Dl Good catch! It's been like that for over six years and you're the first to ever notice. > The odd thing is, if I replace this expression with the one I suggested, the > first set of sweeps (which isn't supposed to move the interface) returns this: > > <eefacjah.png> > The interface then proceeds normally to the right for the rest of the run. I > tried reducing the timestep from 100 to 10 and then 1 and was surprised to > find it became completely unstable and produced a series of overflow errors > (I thought smaller timesteps should always be OK if larger ones converge?). > It was increasing the timestep to 1e4 or larger that finally worked as > intended without moving the interface: > <dgdijibi.png> There are a couple of things going on here, some numerical and some physical, all of which makes this example quite a bit more interesting than as written. First, the numerical: What you're seeing isn't instability, per se, but rather inaccuracy due to exceeding the CFL limit (see http://en.wikipedia.org/wiki/Courant–Friedrichs–Lewy_condition, although that page incorrectly implies that implicit methods completely solve this problem). In short, nothing should convect more than one grid space in one timestep. We can get an estimate of the maximum phase transformation velocity by assuming that the phase changes from 1 to 0 in a single grid spacing and that the diffusivity is Dl at the interface. The term due to the difference in barrier heights is negligible, so: \vec{u}_\phi &\approx \frac{D \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{D \frac{1}{2} V_m}{R T} \left(L_B + L_A) \frac{T_M^A - T_M^B}{T_M^A + T_M^B} \frac{1}{\Delta x} which gives ~0.28 cm/s (we can print max(phaseTransformationVelocity.mag) and see the actual value is 0.14 cm/s, so our estimate is pretty good). In 100 s, that means the concentration wave will try to propagate 28 cm and our simulation domain is only 20 um long, so this is why the solver often can't find a solution. Both 100 s and 10 s are far too big; 100 s happens to converge and 10 s doesn't (more on that later), but both are "wrong". To get a CFL = u * dt / dx < 1, we need a time step of about 1e-5 s. It turns out that we can cheat a little bit, though. Because the phase interface is already in the correct position, we can get a steady-state solution by taking a huge time step as the problem is then almost entirely diffusive. With the corrected diffusivity, I find that dt = 1e4 does the trick. Intermediate time steps are not unstable in the way that explicit diffusion is, but they are inaccurate and they converge poorly, if at all. This is only a guess, but 100 s is more like steady state and 10 s is more like these very inaccurate transients. The steady-state "trick" doesn't work for the second phase of the simulation, where we quench the temperature. Both the phase and concentration profiles are far away from their steady state solutions and so there are large convective driving forces. We can improve the situation a little bit by choosing a different weighting for the diffusivity; if we average the activation energies instead of the diffusivities, then D = Ds**phi * Dl**(1-phi) which tends to make the interface kinetics more like the solid than like the liquid, which reduces the CFL by about an order of magnitude, allowing a corresponding increase in time step. Whether this is physically valid or not is debatable; surface diffusion needn't be like either one, but I would tend to think it was more like liquid state than solid state diffusion. Second, the physical: If you reduce dt to 1e-5 s and run for many time steps, the interface does move to the right. Why is that, when we calculated that the interface should remain stationary? The *equilibrium* solution has the interface in the middle, but that says nothing about the transient path from the initial condition to get to equilibrium. If the diffusivity was uniform, then the concentration wave would propagate at equal rates into solid and liquid and the interface wouldn't move. In this case, the concentration can only change very slowly in the solid compared to the liquid, so the system tends to solidify, concentrating B in the remaining liquid and segregating a lower B concentration in the newly formed solid, consistent with the phase diagram. The higher initial B concentration in the original solid is "frozen in" at this relatively short time scale. It will take about (L/2)**2 / Ds or 1000 s for this excess concentration to dissipate and for the interface to melt back to the middle! ; it just so happens that this is what I empirically found was the time step that gets us close to "steady state". I will need to update the example to reflect these issues. Thanks for bringing this up. _______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
