Jon --
Thanks for this reply. There is something in it I don't understand. Why
is it necessary in fipy constrain the tangential derivative to the solution
to the Poisson equation in fipy? In analytical solutions of Poisson's
equation one can solve the problem either by specifying the normal
derivative of the solution on the boundary OR the value of the function on
the boundary. To do both generally over-specifies the solution.
Specifying the tangential gradients is equivalent to specifying the
value of the solution on the boundary (to within a constant), for one can
integrate the tangential gradient along the boundary to find the solution
along the boundary.
Why is it in general consistent to specify both the normal and
tangential boundary conditions?
Thanks,
Jamie
On Fri, May 5, 2017 at 3:33 PM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:
> Jamie -
>
>
>
> Sorry it took so long to get back to this. I'd actually figured out what
> the issue was several days ago, but I haven't had time to verify and write
> anything up until now.
>
>
>
> The short answer is that you need to constrain the gradient of PsiInvert
> to be equal to the gradient of PsiTruth, but what you're doing (aside from
> being a lot of unnecessary typing and I'm pretty sure having low-order
> spacial accuracy) is only constraining the normal component of the boundary
> gradient. You need the tangential component, too.
>
>
>
> I've updated your script below with the proper boundary gradients, e.g.,
>
>
>
> psiInvert.faceGrad.constrain(psiTrue.faceGrad,where=yeq0bound)
>
>
>
> and the recommended way to set omega to the Laplacian:
>
>
>
> omega = psiTrue.faceGrad.divergence
>
>
>
> I get errors of O(1e-10) for all three cases.
>
>
>
> In the course of working through this, I realized that I'd not finished
> converting .faceGrad to math in what I posted before, such that I also
> neglected the tangential term. I've updated the notebook at
>
>
>
> https://urldefense.proofpoint.com/v2/url?u=https-3A__gist.
> github.com_guyer_6b7699531f75b353f49a0cf36d683a
> a4&d=DwIGaQ&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=
> mjWx1u6TsXKmtow-qItsKS4VE7r1nIMW-c395mi9wzM&s=EQ1vHADeK-
> 0hEbWdEYmXbmGl2msCVM8sAPo2PX5wb7w&e=
>
>
>
> to correct this and also to illustrate a *very hacky* implementation of
> boundary upwinding.
>
>
>
> - Jon
>
>
>
> > On Apr 26, 2017, at 9:02 AM, James Pringle <jprin...@unh.edu> wrote:
>
> >
>
> > Dear all --
>
> >
>
> > I think I have discovered a bug in the implementation of Neumann
> boundary conditions in fipy for the solution of Poisson's equation. If it
> is not a bug, it exposes a significant inefficiency in the numerics.
> Attached is a test code (DemoBug.py) which illustrates the problem. The
> test code takes a given function PsiTruth, takes the the Laplacian of it,
> and then solves Possion's equation to recover Psi (we call this recovered
> version of Psi "PsiInvert," and it should ideally be exactly equal to
> PsiTruth):
>
> >
>
> > \Del^2 PsiInvert = \Del^2 PsiTruth
>
> >
>
> > The problem is solved in a rectangular domain of size (Lx,Ly) with the
> following boundary condition's:
>
> > • PsiInvert=PsiTruth(x,Ly) on y=Ly (a Dirchelet BC)
>
> > • The normal derivative of PsiInvert = The normal derivative of
> PsiTruth on all other boundaries.
>
> > The normal gradient of PsiTruth for the boundary condition can be solved
> either directly from PsiTruth or by using the fipy.faceGradAverage() on
> PsiTruth on the mesh. This choice is made in the code after the comment
> "#Choose Boundary Condition". I have three choices of the PsiTruth to
> examine. The error is defined as the STD(PsiInvert-PsiTruth). All have a
> peak amplitude of about 1.0, so the standard deviation of the error in the
> solution should be <<1.
>
> > • An isolated gaussian that does not extend to the boundaries;
> this works very well with a STD(error) of 3.7e-4
>
> > • cos(2*pi/wavelength*x)*cos(2*pi/wavelength*y), where wavelength
> is 160 the resolution of the grid (so is very well resolved). This works
> less well but ok with an error of 5.0e-2 (two orders of magnitude larger
> than above, about 5% of the solution magnitude).
>
> > • abs(cos(2*pi/wavelength*x)*cos(2*pi/wavelength*y)), same wave
> length as above. Error is 4.8e-1; the error in the solution is 3 orders of
> magnitude greater than case 1 and an order of magnitude worse than case 2.
> It is 50% of the solution magnitude.
>
> > In all of the results above, the boundary normal gradients for the BC
> were computed using faceGradAverage() on PsiTruth. Worryingly, the errors
> in case 2 and 3 are an order of magnitude greater if the true gradient is
> computed directly from the underlying function for PsiTruth!
>
> >
>
> > If you look at the solution (PsiInvert_Case3.png) and and error
> (Error_Case3.png) to case 3, attached as images, you will see that the
> error is a large scale feature, essentially a plane tilting upward with
> increasing x over the whole domain. In the similar plots for Case 2 you can
> see that the magnitude of the error increases and decreases as the
> underlying function becomes positive and negative. The net results of these
> fluctuations is that the magnitude of the error decreases about 1
> wavelength into the interior, as one would expect for an oscillating
> boundary condition applied to a Laplacian.
>
> >
>
> > A hint as to the origin of the error can be found by looking at the
> normal derivative at the boundary x=Lx
> (BoundaryNormalDerivative_Case[2,3].png).
> I plot the normal and tangential derivative directly calculated from the
> function that defines PsiTrue to the gradients calculated by
> faceGradAverage(). You can see that there is an error in the normal
> derivative calculated by faceGradAverage of approximately 2; however even
> if the correct normal gradient is used as the boundary condition, the error
> is even larger! This suggests their is something deeply suboptimal in the
> implementation of the boundary condition.
>
> >
>
> > Now, the error does go towards zero as the resolution is increased; The
> error in case 3 is reduced a factor of 8 to about 6% when the resolution is
> increased by a factor of 8; in general the error scales linearly with
> resolution. But this means that to get acceptable performance (about a 5%
> error), the resolution must be 1280 smaller than the wavelength of the
> problem to be solved. This is extreme, and suggests something deeply
> non-optimal.
>
> >
>
> > Roughly, the error scales as the (resolution)*(gradient at boundary); if
> the normal gradient is 0, there is little error.
>
> >
>
> > The same errors exist if the domain is decomposed into triangles instead
> of quadrilaterals. Please let me know if you would like any other test
> cases or work on my part.
>
> >
>
> > I would very much appreciate any help into how to eliminate or reduce
> this error. It is much more than I expect from similar solver codes I have
> used in Fortran (admittedly, finite difference codes like mud2).
>
> >
>
> > Thank you very much,
>
> > 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
>
> >
>
> >
>
> >
>
> >
>
> >
>
> > <BoundaryNormalDerivative_Case2.png><BoundaryNormalDerivative_
> Case3.png><Error_Case2.png><Error_Case3.png><PsiInvert_
> Case2.png><PsiInvert_Case3.png><TruePsi_Case2.png><
> TruePsi_Case3.png><DemoBug.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=mjWx1u6TsXKmtow-qItsKS4VE7r1nIMW-c395mi9wzM&s=
> PIUFyqWI6sR6wZweLlc5veaRvK3MuQ_UrCfVJ88iUA4&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=
> mjWx1u6TsXKmtow-qItsKS4VE7r1nIMW-c395mi9wzM&s=qNyS7ky4clCpAGEgwm2pF5lXe-
> LHjqe1hph5TZLaFd4&e= ]
>
>
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
