Hi Jonathan, Thanks for the reply, and you're right. Most of the simulations I'd been doing involved only positive quantities like concentration, so I was safe using harmonicFaceValue, but my displacement variables u1 and u2 can be positive or negative. I just re-ran it using divergence operators on arithmetic face values of displacements to confirm, and with that it works fine.
Thanks, Ray On Tue, Jun 25, 2013 at 12:16 PM, Jonathan Guyer <[email protected]> wrote: > > On Jun 24, 2013, at 4:03 PM, Raymond Smith <[email protected]> wrote: > > I'm having confusion about what the divergence operator is doing in my > situation. I ran two simulations and got very different results using two > similar scripts, attached. The scripts are relatively long, but I'm > including the full scripts in case there's anything else that someone may > end up wanting to reference. My focus now is on the "diff" of the two > scripts (lines 145-174 and lines 196-225). Basically, I'm not sure why > they're different. I'm in 2 dimensions with a cellVariable called c, and > the term I want in my equation is > elem1 * dc/dx + elem2 * dc/dy, > where elem1 and elem2 are known, constant parameters. > > > My guess would have been that the treatments should be equivalent (except > perhaps in the orders of accuracy?), but the results (movie_DIV.mpg for > *divop.py and movie_DOT.mpg for the *dotop.py file) are very dissimilar. > > > Is there something I'm missing about how the divergence operator is acting > on the FaceVariables? > > > It appears to me that the difference is not DIV vs. DOT, per se, but > rather your use of the harmonicFaceValue for fields that are not uniformly > positive. The harmonic average is only meaningful for quantities that are > greater than zero. > > > >>> import fipy as fp > >>> from fipy import numerix as nx > >>> mesh = fp.Grid2D(nx=100, Lx=4 * nx.pi, ny=100, Ly=4 * nx.pi) > >>> x, y = mesh.cellCenters[0], mesh.cellCenters[1] > >>> c = y * nx.sin(y) * nx.cos(x) > >>> fp.Viewer(c).plot(filename="field.png") > > >>> elem1 = 3. > >>> elem2 = 5. > > Note, there's no reason to make the nhats and dot them. Just ask for the > scalar component you want: > > >>> dotty = elem1 * c.grad[0] + elem2 * c.grad[1] > >>> fp.Viewer(dotty).plot(filename="dotty.png") > > > >>> a1 = fp.Variable((elem1, elem2)) > > The harmonic mean produces some surprising artifacts: > > >>> divvy_harmonic = (a1 * c.harmonicFaceValue).divergence > >>> fp.Viewer(divvy_harmonic).plot(filename="divvy_harmonic.png") > > The arithmetic mean does not: > > >>> divvy_arithmetic = (a1 * c.arithmeticFaceValue).divergence > >>> fp.Viewer(divvy_arithmetic).plot(filename="divvy_arithmetic.png") > > The difference between the divergence of the arithmetic mean and the sum > of the components is quite small: > > >>> fp.Viewer(divvy_arithmetic - dotty).plot(filename="diff.png") > > > Also, out of curiosity, examples.diffusion.anisotropy uses > fp.numerix.NUMERIX.dot for some dotting. Is there a reason for the > preference of that over fp.numerix.dot, and if so, how is it different? > > > fp.numerix.dot performs dot products on each position of a field of > vectors. > fp.numerix.NUMERIX.dot (an alias for numpy.dot) does not understand how > FiPy treats fields of vectors, but tries to dot the entire field as a > single vector (or tensor). > > This is why we always stress to use the functions from fipy.numerix, > rather than from numpy directly. > > numpy.dot is, nonetheless, a useful function for FiPy's internals. > > > > > _______________________________________________ > 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 ]
