I haven't had a chance to look at your second message. I'll get back to you when I've gone through it.
> On Apr 26, 2017, at 2:11 PM, James Pringle <[email protected]> wrote: > > 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) > <[email protected]>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_6b7699531f75b353f49a0cf36d683aa4&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 <[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] > > 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 > [email protected] > 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 > [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 ]
