Jon -- Thank you for that clear explanation. What remains unclear to me is why the .faceGrad.constrain() is not constraining the solution to Poisson's equation properly. I have given a clearer example of that in a follow on email titled "bug in Neumann boundary condition to Laplacian"
Thanks, Jamie On Wed, Apr 26, 2017 at 12:06 PM, Guyer, Jonathan E. Dr. (Fed) < jonathan.gu...@nist.gov> wrote: > 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://urldefense.proofpoint.com/v2/url?u=https-3A__gist. > github.com_guyer_6b7699531f75b353f49a0cf36d683a > a4&d=DwIGaQ&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m= > 8Otm5er2tke4qe_pfoFqjD7xv-xYoMTBdUD7HW01qxg&s=mIlZMfmcm_ > vIrOBVFcGZFNwOGxucRHOzR7wyx2D946U&e= > > 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 <jprin...@unh.edu> 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 <james.prin...@unh.edu> > 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 > > fipy@nist.gov > > https://urldefense.proofpoint.com/v2/url?u=http-3A__www. > ctcms.nist.gov_fipy&d=DwIGaQ&c=c6MrceVCY5m5A_KAUkrdoA&r= > 7HJI3EH6RqKWf16fbYIYKw&m=8Otm5er2tke4qe_pfoFqjD7xv-xYoMTBdUD7HW01qxg&s= > tWYZmHLXWb-PM9h_Df4nqHKzLVCScRwAaVjiV5AXVI0&e= > > [ NIST internal ONLY: https://urldefense.proofpoint. > com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_ > fipy&d=DwIGaQ&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m= > 8Otm5er2tke4qe_pfoFqjD7xv-xYoMTBdUD7HW01qxg&s=McRS8G7i- > 4lzkmHEnKM74LmRScSRxZ11SzYObpGmTXg&e= ] > > > _______________________________________________ > fipy mailing list > fipy@nist.gov > https://urldefense.proofpoint.com/v2/url?u=http-3A__www. > ctcms.nist.gov_fipy&d=DwIGaQ&c=c6MrceVCY5m5A_KAUkrdoA&r= > 7HJI3EH6RqKWf16fbYIYKw&m=8Otm5er2tke4qe_pfoFqjD7xv-xYoMTBdUD7HW01qxg&s= > tWYZmHLXWb-PM9h_Df4nqHKzLVCScRwAaVjiV5AXVI0&e= > [ NIST internal ONLY: https://urldefense.proofpoint. > com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_ > fipy&d=DwIGaQ&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m= > 8Otm5er2tke4qe_pfoFqjD7xv-xYoMTBdUD7HW01qxg&s=McRS8G7i- > 4lzkmHEnKM74LmRScSRxZ11SzYObpGmTXg&e= ] >
_______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]