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://gist.github.com/guyer/6b7699531f75b353f49a0cf36d683aa4
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
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
from pylab import *
from numpy import *
import time
import matplotlib.tri as tri
import fipy
#this code creates and arbitrary streamfunction psi, takes the
#Laplacian of it, and then inverts the Laplacian back by solving
#Poisson's equation. The answer should be thw same as the test
#streamfunction.
#==================================================================
#define domain with Grid2D; it is rectangular, with dimension Lx by Ly
dx=1.0e0
dy=dx
Ly=500e0
Lx=110.0e0;
nx=int(round(Lx/dx)); ny=int(round(Ly/dy))
mesh=fipy.Grid2D(dx=dx,dy=dx,nx=nx,ny=ny)
assert max(mesh.faceCenters.value[0,:])==Lx,'Lx must be multiple of dx, or the BCs will fail'
assert max(mesh.faceCenters.value[1,:])==Ly,'Ly must be multiple of dy, or the BCs will fail'
#==================================================================
#define the stream function we use to test. The first is solved
#perfectly. The solution to thesecond has a small error, and the
#third case leads to a very inaccurate solution.
#Choose Case
Case=3
if Case==1 : #Case 1, an isolated gaussian that does not touch borders. This case works perfectly.
print 'Running Case 1'
rPsi=12.0
Psi0=1.0
def psiFunc(x,y):
return Psi0*exp(-(x-Lx/2)**2/(2*rPsi**2))*exp(-(y-Ly/2)**2/(2*rPsi**2))
elif Case==2: # Case 2, This case sort of works;
print 'Running Case 2'
wavelength=80.0e0*2
print 'wavelength/dx=',wavelength/dx
def psiFunc(x,y):
return cos(2*pi/wavelength*x)*cos(2*pi/wavelength*y)
elif Case==3: #Case 3, this leads to massive errors
print 'Running Case 3'
wavelength=80.0e0*2
print 'wavelength/dx=',wavelength/dx
def psiFunc(x,y):
return abs(cos(2*pi/wavelength*x)*cos(2*pi/wavelength*y))
else:
assert False,'you got to pick a case that exists!'
psiTrue=fipy.CellVariable(name='PsiTrue',mesh=mesh,
value=psiFunc(mesh.x,mesh.y))
#draw true solution
figure(1)
clf()
tricontourf(mesh.x,mesh.y,psiTrue.value,30)
colorbar()
title('True Psi')
savefig('TruePsi.png')
#==================================================================
#Calculate Laplacian of Psi
#start by calculating divergence of stream function
print 'Calculating Laplacian in interior with faceGrad.divergence'
omega = psiTrue.faceGrad.divergence
#==================================================================
#Now invert omega for Psi using \Del^2 Psi=omega.
#First set up boundary conditions
#The boundary conditions are from the velocity on the faces
#
#Psi_x=Grad_x on the y=0 boundaries and Psi_y=Grad_y on the x=0,Lx boundary.
#
#At y=Ly, set Psi BC to Psi(at y=Ly)
#
#first, find locations of boundaries
xFace,yFace=mesh.faceCenters
yeq0bound=yFace==0.0
yeq1bound=yFace==Ly
xeq0bound=xFace==0.0
xeq1bound=xFace==Lx
#because these boundaries are specified on the faces, we need to
#recalculate Grad_y and Grad_x on the faces
#calculate gradient on face from fipy
#'Calculating gradient of psi on faces for BC with faceGradAverage'
jnk=psiTrue.faceGradAverage;
Grad_xface_fipy=jnk[0,:]
Grad_yface_fipy=jnk[1,:]
#calculating gradient of psi on faces quasi-analytically
Dx=dx/2
Grad_xface_qa=(psiFunc(xFace+0.5*Dx,yFace)-psiFunc(xFace-0.5*Dx,yFace))/Dx
Grad_yface_qa=(psiFunc(xFace,yFace+0.5*Dx)-psiFunc(xFace,yFace-0.5*Dx))/Dx
#compare the two gradient calculations on a boundary
#where=isfinite(xFace); titleStr='all'
where=xeq1bound; titleStr='x=Lx'
figure(4)
clf()
subplot(2,1,1)
plot(Grad_xface_qa[where],Grad_xface_fipy[where],'r.',alpha=0.2,mec='none',label='true vrs fipy');
plot(Grad_xface_qa[where],Grad_xface_qa[where],'k-',mec='none',label='line of y=x');
legend()
ylabel(r'$\partial \Psi/\partial x$ from fipy'),xlabel(r'$\partial \Psi/\partial x$');grid()
title('Gradient on '+titleStr+' correct gradient vrs. that calculated by fipy')
subplot(2,1,2)
plot(Grad_yface_qa[where],Grad_yface_fipy[where],'r.',alpha=0.2,mec='none');
plot(Grad_yface_qa[where],Grad_yface_qa[where],'k-',mec='none');
ylabel(r'$\partial \Psi/\partial y$ from fipy'),xlabel(r'$\partial \Psi/\partial y$');grid()
savefig('BoundaryNormalDerivative.png')
#Choose Boundary Condition
#use boundary condition as calculated by fipy, or by manually differentiating PsiTruth
if True:
print 'Calculating gradient of psi on faces for BC with faceGradAverage'
Grad_xface=Grad_xface_fipy
Grad_yface=Grad_yface_fipy
else:
print 'Calculating gradient of psi on face for BC quasi-analytically'
Grad_xface=Grad_xface_qa
Grad_yface=Grad_yface_qa
#define boundary conditions and solve
psiInvert=fipy.CellVariable(name='psiInvert',mesh=mesh,value=0.0)
faceNorm=mesh.faceNormals
psiInvert.faceGrad.constrain(psiTrue.faceGrad,where=yeq0bound)
psiInvert.constrain(psiTrue.faceValue,where=yeq1bound)
psiInvert.faceGrad.constrain(psiTrue.faceGrad,where=xeq0bound)
psiInvert.faceGrad.constrain(psiTrue.faceGrad,where=xeq1bound)
#=====================================================================
#set up differential equation and solve
#now define variable and set BCs
#solve \Del^2 Psi=omega
tic=time.time()
sourceTermPsi=fipy.CellVariable(name='psiForcing',mesh=mesh,value=-omega)
DiffCoefPsi=fipy.CellVariable(mesh=mesh,rank=0,value=1.0)
eqPsi=(fipy.DiffusionTerm(var=psiInvert,coeff=DiffCoefPsi)+sourceTermPsi)
eqPsi.solve(var=psiInvert)
#=======================================================================
#now print statistics
print 'Done finding psi in',time.time()-tic
print 'max and min solved for Psi are %e,%e'%(amin(psiInvert.numericValue),amax(psiInvert.numericValue))
print 'STD of error between true and inverted for psi is %e'%(std((psiInvert.value-psiTrue.value)),)
print 'STD of PsiTrue is %f, solved for psi %f'%(std(psiTrue.value),std(psiInvert.value))
#=======================================================================
#now plot error in solution and solution
figure(2)
clf()
tricontourf(mesh.x,mesh.y,(psiInvert.value)-(psiTrue.value),30)
colorbar()
title('Error=(Solved for Psi)-(True Psi)')
savefig('Error.png')
figure(3)
clf()
tricontourf(mesh.x,mesh.y,psiInvert.value,30)
colorbar()
title('Psi found by solving Poissons eq.')
savefig('PsiInvert.png')
draw()
show()
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
