James - Thank you for your test script. It helped me understand what you were seeing and what behavior you were looking for. I've put together a notebook that I think concisely illustrates the source of the behavior you're encountering.
https://gist.github.com/guyer/6b7699531f75b353f49a0cf36d683aa4 The bottom line is that FiPy assumes zero gradient on all boundaries unless some constraint is applied. FiPy does not do any upwinding to project interior values and slopes out to the boundary. Further, Dirichlet constraints are reflected in boundary gradients, but Neumann constraints are not reflected in boundary values. Doing the latter (while desirable) would require unwinding (expensive) and care to avoid circular definitions (expensive). We'll think about how to achieve this and we'd certainly welcome code contributions from anybody who wants to tackle it. Finally, take care with your parentheses. .grad, .faceGrad, etc. are properties, not methods. Appending () to the end causes them to be evaluated to static numeric values (phi.grad() is the same as phi.grad.value). - Jon > On Apr 21, 2017, at 6:33 PM, James Pringle <[email protected]> wrote: > > To clarify and answer my question, and to give an analysis code for those > interested in figuring these out, here is some text. (I would love comments > from those who know more!) > > First a restatement of the question. To test the solvers, create a field > PsiTruth, and solve for PsiInvert with > > \Del^2 PsiInvert = \Del^2 PsiTruth > > With various boundary conditions. The attached code does this. (To understand > the notation in the code, see my prior email). My conclusions are that for > maximum accuracy: > • Compute \Del^2 PsiTruth in the interior on center points with the > grad() operator. It is significantly more consistent with the solver numerics > than the leastSquaresGrad() operator. This results in smaller error. > • On the boundaries, use .faceGradAverage() to calculate the gradients > of PsiTruth; .faceGrad() is incorrect on the boundaries. > • At least part of one boundary must be of the form > PsiInvert.constrain(), for otherwise the solution is arbitrary to within an > additive constant, and that constant can be so large that it compromises the > accuracy of the solution due to the finite number of significant digits in > floating point numbers. > • When gradient boundary conditions are specified, accuracy increases > linearly with 1/resolution, and is much less than when the boundary values > are fixed. > I hope this helps someone, and I would like to hear from others about the > choice of gradient operators and how to achieve greater accuracy in > calculating and specifying gradients on the boundary. > > In particular, I would love to have a gradient calculation that is exactly > consistent with how the gradient is specified in the boundary condition > constraint. > > It would also be nice to have a Laplacian operator that is consistent with > the solver in the model -- or do I get that by using the .grad() operator > twice? > > Thanks all, > Jamie > > > On Fri, Apr 21, 2017 at 11:28 AM, Pringle, James <[email protected]> > wrote: > Dear all -- > > I am using fipy to solve a series of fluid dynamics problems. As part of > this, I need to compute a stream function Psi from the vorticity of flow > omega=V_x-U_y. U and V are calculated from the gradient of the solution to > another set of elliptic equations, and are available on the same mesh as Psi > will be calculated on. > > The equation relating the stream function Psi to the vorticity omega is > > \Del^2 Psi = omega > > And the boundary conditions on Psi are that the gradient normal to the > boundary is equal to the flow parallel to the boundary. All very standard. > > I calculate velocities U and V from the gradient of another variable (a > pressure like term). My question is this: which of the routines that > calculate the gradient of a fipy variable on the faces is most consistent > with how the normal derivative is specified on the boundary? Which of > .faceGrad() or faceGradAverage() is most consistent with how the gradient is > implemented in Psi.faceGrad.constrain()? or some other gradient calculation? > I ask because .faceGrad() does not seem to calculate the gradient on the > boundary faces... > > Thanks, > Jamie Pringle > > Director Ocean Process Analysis Laboratory in the Insitute for Earth, Ocean > and Space > Associate Professor of Earth Sciences > University of New Hampshire > http://oxbow.sr.unh.edu > > <Psi_inversion_test.py>_______________________________________________ > 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 ]
