I don't know. The form you write down in "split plan 1" is what we write as eq3, but you can see that we don't use eq3 in our solution. As I said the other day, the form of eq2 in our example is chosen because that Taylor expansion has generally better convergence properties than other forms of that source.
There are a multitude of ways to decompose the problem. Your "split plan 2" might be better. Run time is one metric, although not a great one. I'd look at the convergence rate on the residuals, as well as compare the solutions from the different schemes to see if they're solving to the same accuracy. On Feb 23, 2015, at 4:47 AM, Ronghai Wu <[email protected]> wrote: > If "psi = epsilon**2*laplace phi"(Split plan 2 in codes) is a better choice > than "psi = dfdphi - epsilon**2*laplace phi"(Split plan 1 in codes)? Because > "Split plan 2" runs much faster than "Split plan 1". Is there any hidden > trouble using "Split plan 2"? > > > import fipy as fp > nx = 64 > ny = 64 > mesh = fp.Grid2D(nx=nx, ny=ny, dx=1., dy=1.) > > psi = fp.CellVariable( hasOld=1,mesh=mesh ) > phi = fp.CellVariable( mesh = mesh, value = > fp.GaussianNoiseVariable(mesh=mesh, mean=0.5, > variance=0.01).value, hasOld=1 ) > > > D = a = 1. > epsilon = 0.5 > dfdphi_ = a**2 * (1 - phi) * (1 - 2 * phi) > d2fdphi2 = a**2 * (1 - 6 * phi * (1 - phi)) > > #Split plan 1 > eq1 = (fp.TransientTerm(coeff=1.,var=phi) == > fp.ImplicitDiffusionTerm(coeff=D, var=psi)) > > eq2 = (fp.ImplicitSourceTerm(coeff=1., var=psi) > == fp.ImplicitSourceTerm(coeff=dfdphi_, var=phi) - > fp.DiffusionTerm(coeff=epsilon**2, var=phi)) > > > #Split plan 2 > #eq1 = (fp.TransientTerm(coeff=1.,var=phi) == > fp.ImplicitDiffusionTerm(coeff=d2fdphi2, var=phi) + > fp.ImplicitDiffusionTerm(coeff=-D, var=psi)) > > #eq2 = (fp.ImplicitSourceTerm(coeff=1., > var=psi) == fp.DiffusionTerm(coeff=epsilon**2, var=phi)) > > > eq = eq1 & eq2 > dt = 0.5 > viewer = fp.Viewer(vars=(phi, psi)) > for t in range(1000): > phi.updateOld() > psi.updateOld() > res = 1.e5 > while res > 1.e-1: > res = eq.sweep(dt=dt) > print res > viewer.plot() > > > > > > > > > > On 20/02/15 19:22, Guyer, Jonathan E. Dr. wrote: >> On Feb 20, 2015, at 6:39 AM, Ronghai Wu <[email protected]> >> wrote: >> >> >>> (1) The 4th-order is split by "psi = d2fdphi2(phi-phiold) + dfdphi - >>> epsilon**2*laplace phi". I do not understand why we need >>> d2fdphi2(phi-phiold)? >>> >> See the discussion on linearization of the source in >> >> >> http://www.ctcms.nist.gov/fipy/examples/phase/generated/examples.phase.simple.html >> >> >> You could write "psi = dfdphi - epsilon**2*laplace phi", but dfdphi would be >> a fully explicit source. By expanding >> >> source = source_old + d(source)/d(phi) (phi - phi_old) >> >> you get a more implicit version of the same term, with better convergence >> properties. >> >> >> _______________________________________________ >> fipy mailing list >> >> [email protected] >> http://www.ctcms.nist.gov/fipy >> >> [ NIST internal ONLY: >> https://email.nist.gov/mailman/listinfo/fipy >> ] >> > > -- > ------------------------------------------ > Ronghai Wu > > Institute of Materials Simulation (WW8) > Department of Materials Science and Engineering > University of Erlangen-Nürnberg > Dr.-Mack-Str. 77, 90762 Fürth, Germany > > Tel. +49 (0)911 65078-65064 > > _______________________________________________ > 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 ]
