Thanks for the explanation, I think I understand now. Adam
On 11/11/2012 12:15 PM, Jonathan Guyer wrote: > 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 midd! le! > ; 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 ] _______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
