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 ]
