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:
1. 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.
2. On the boundaries, use .faceGradAverage() to calculate the gradients
of PsiTruth; .faceGrad() is incorrect on the boundaries.
3. 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.
4. 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
>
from pylab import *
from numpy import *
import time
import matplotlib.tri as tri
import fipy
#this code creates and arbitrary streamfunction psi, finds velocity
#from it, and then inverts the velocity back to the stream function
#using the vorticity equation. The result should be the same as the
#original streamfunction, to within an additive constant.
#==================================================================
#define domain
dx=1.0
Lx=125.0
nx=int(round(Lx/dx)); ny=int(round(Lx/dx))
mesh=fipy.Grid2D(dx=dx,dy=dx,nx=nx,ny=ny)
assert max(mesh.faceCenters.value[0,:])==Lx,'dx must be multiple of Lx, or the BCs will fail'
#define streamfunction as cos(2*pi/wavelength*mesh.x)*cos(2*pi/wavelength*mesh.y) and plot
wavelength=80.0
psiTrue=fipy.CellVariable(name='PsiTrue',mesh=mesh,
value=cos(2*pi/wavelength*mesh.x)*cos(2*pi/wavelength*mesh.y))
print 'wavelength/dx=',wavelength/dx
figure(1)
clf()
tricontourf(mesh.x,mesh.y,psiTrue.value,30)
colorbar()
title('True Psi')
#==================================================================
#calculate U and V from streamfunction with U=-Psi_y and V=Psi_x,
#and plot with quiver.
if False:
#calculate gradient with leastSquaresGrad
print 'Calculating gradient in interior with leastSquaresGrad'
jnk=psiTrue.leastSquaresGrad()
Vcenter=jnk[0,:]
Ucenter=-jnk[1,:]
#==================================================================
#Now calculate the vorticity omega=V_x-U_y
Ufipy=fipy.CellVariable(name='U',mesh=mesh,value=Ucenter); Ugrad=Ufipy.leastSquaresGrad()
Vfipy=fipy.CellVariable(name='V',mesh=mesh,value=Vcenter); Vgrad=Vfipy.leastSquaresGrad()
omega=Vgrad[0,:]-Ugrad[1,:] #V_x-U_y
elif True:
#calculate gradient with Grad
print 'Calculating gradient in interior with Grad'
jnk=psiTrue.grad()
Vcenter=jnk[0,:]
Ucenter=-jnk[1,:]
#==================================================================
#Now calculate the vorticity omega=V_x-U_y
Ufipy=fipy.CellVariable(name='U',mesh=mesh,value=Ucenter); Ugrad=Ufipy.grad()
Vfipy=fipy.CellVariable(name='V',mesh=mesh,value=Vcenter); Vgrad=Vfipy.grad()
omega=Vgrad[0,:]-Ugrad[1,:] #V_x-U_y
else:
#analytical solution
print 'Calculating gradient in interior analytically'
omega=-2*((2*pi/wavelength)**2.0)*cos(2*pi/wavelength*mesh.x)*cos(2*pi/wavelength*mesh.y)
if False: #plot velocity
nxcull=15; quiver(mesh.x[::nxcull],mesh.y[::nxcull],
Ucenter[::nxcull],Vcenter[::nxcull],scale=2.0e0)
#==================================================================
#Now invert omega for Psi using \Del^2 Psi=omega.
#The boundary conditions are from the velocity on the faces
#
#Psi_x=V on the y=0,Lx boundaries and Psi_y=-U on the x=0,Lx boundary.
#
#because these boundaries are specified on the faces, we need to
#recalculate U and V on the faces
if True:
print 'Calculating gradient of psi on faces for BC with faceGradAverage'
jnk=psiTrue.faceGradAverage();
Vface=jnk[0,:]
Uface=-jnk[1,:]
elif True:
#NOTE WELL -- this is a very bad solution, since it gives a normal derivative of 0
#regardless of psiTrue
print 'Calculating gradient of psi on faces for BC with faceGrad -- which is wrong'
jnk=psiTrue.faceGrad();
Vface=jnk[0,:]
Uface=-jnk[1,:]
else:
#instead use analytical values for the boundary condition,
#again with U=-Psi_y, V=Psi_x
print 'Calculating gradient in psi on faces for BC with analytical solutions'
xFace,yFace=mesh.faceCenters
Vface=-2*pi/wavelength*sin(2*pi*xFace/wavelength)*cos(2*pi*yFace/wavelength)
Uface= 2*pi/wavelength*cos(2*pi*xFace/wavelength)*sin(2*pi*yFace/wavelength)
#nxcull=15; quiver(xFace[::nxcull],yFace[::nxcull],
# Uface[::nxcull],Vface[::nxcull],scale=2.0e0,color='r')
#set up differential equation and solve
#first, find locations of boundaries
xFace,yFace=mesh.faceCenters
yeq0bound=yFace==0.0
yeq1bound=yFace==Lx
xeq0bound=xFace==0.0
xeq1bound=xFace==Lx
#at some location, we need to constrain psi, and not just the normal
#gradient of psi, at the boundary. Otherwise it will differ by an
#additive constant. Boundary points where the value is constrained are
#defined by the vector pntsTofix
pntsToFix=copy(xeq1bound) #just fix on the x=Lx boundary
#pntsToFix=copy(yeq0bound) #just fix on the y=0 boundary
#pntsToFix=xeq1bound | xeq0bound | yeq1bound | yeq0bound #fix all boundaries
#remove fixed points from points which are fixed by gradient
xeq0bound[pntsToFix]=False; xeq1bound[pntsToFix]=False
yeq0bound[pntsToFix]=False; yeq1bound[pntsToFix]=False
fixValueHere=isnan(yFace) # all false
fixValueHere[pntsToFix]=True
#now define variable and set BCs
psiInvert=fipy.CellVariable(name='psiInvert',mesh=mesh,value=0.0)
faceNorm=mesh.faceNormals
psiInvert.faceGrad.constrain(-Uface*faceNorm,where=yeq0bound)
psiInvert.faceGrad.constrain(-Uface*faceNorm,where=yeq1bound)
psiInvert.faceGrad.constrain(Vface*faceNorm,where=xeq0bound)
psiInvert.faceGrad.constrain(Vface*faceNorm,where=xeq1bound)
#set value at some location on boundary, so constant offset is constrained
psiInvert.constrain(psiTrue.faceValue,fixValueHere)
#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,solver=fipy.LinearLUSolver())
#eqPsi.solve(var=psiInvert)
print 'Done finding psi in',time.time()-tic
print 'max and min psi are %e,%e'%(amin(psiInvert.numericValue),amax(psiInvert.numericValue))
#print 'Psi at fixed point in true %f and inverted %f Psi',valueAtOnePoint,psiInvert.faceValue.value[fixValueHere][0]
#when comparing, it is important to remove means, because the solution
#is undefined to within a constant
meanPsiTrue=mean(psiTrue.value)
meanPsiInvert=mean(psiInvert.value)
print 'STD of error is %e'%(std(((psiInvert.value-meanPsiInvert)-(psiTrue.value-meanPsiTrue))),)
print 'STD of PsiTrue is %f, PsiInvert %f'%(std(psiTrue.value),std(psiInvert.value))
figure(2)
clf()
tricontourf(mesh.x,mesh.y,(psiInvert.value-meanPsiInvert)-(psiTrue.value-meanPsiTrue),30)
colorbar()
title('(Inverted Psi)-(True Psi), both with means removed')
figure(3)
clf()
tricontourf(mesh.x,mesh.y,psiInvert.value,30)
colorbar()
title('Inverted Psi')
draw()
show()
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
