Thank you, Jonathan, for the insight. Sorry about my confusion -- I'm still
slowly teaching myself the fundamentals as a biologist.
I've now read more on the topic. Based on your comments:
"
> By constraining var.faceGrad to zero, you are saying that the boundary
> flux is
> \phi b \tanh(alpha*r) (pos-den)/r
> "
>
This is not a conservative boundary condition.
it makes more sense why having *var.faceGrad.constrain(0,
where=m.exteriorFaces)* causes the solution to blow up if advection points
rightward toward a point-attractor ("den") and to leak if advection points
leftward.
I think I've now also correctly set up the Robin BC (line 38 in attached
script) -- the solution's integral over space ("cell volume") is indeed
conserved at 1 throughout the simulation. However, not defining my boundary
condition at all appears to give exactly the same results. Is this because
FiPy assumes Robin BC if given no specification?
Thanks,
Yun
On Thu, Aug 25, 2016 at 10:52 AM, Guyer, Jonathan E. Dr. (Fed) <
[email protected]> wrote:
> I'm not sure why you conclude from that example that
>
> var.faceGrad.constrain(0, where=m.exteriorFaces)
>
> only applies to the advection component. Where does it say that?
>
>
> In the code you sent, with
>
> eq3 = TransientTerm() == DiffusionTerm(coeff=D)+
> convection(coeff=faceVelocity)
>
> the flux is
>
> D \nabla \phi + \phi b \tanh(alpha*r) (pos-den)/r
>
> By constraining var.faceGrad to zero, you are saying that the boundary
> flux is
>
> \phi b \tanh(alpha*r) (pos-den)/r
>
> This is not a conservative boundary condition.
>
> > On Aug 24, 2016, at 8:59 PM, Yun Tao <[email protected]> wrote:
> >
> > Hi Dan,
> >
> > Wow. I did not expect that at all, thanks for that important
> information! I just tried searching for ways to implement the Robin BCs in
> FiPy. This site example <http://www.ctcms.nist.gov/fip
> y/examples/convection/generated/examples.convection.robin.html> appears
> to show that the var.faceGrad.constrain(0, where=m.exteriorFaces) command I
> mentioned earlier handles only the advection component, while, I presume,
> the diffusion component is subject to Neumann BCs by default. If that is
> so, then it's even more puzzling why having that line in my code blows up
> the solution, yet the solution seems well-behaved under just the Neumann
> BCs. Am I missing something else here?
> >
> > Thanks,
> > Yun
> >
> > On Wed, Aug 24, 2016 at 8:24 PM, Daniel Farrell <[email protected]>
> wrote:
> > Hello Yun,
> >
> > I just briefly looked at the code. Seems like you are solving something
> like an advection diffusion problem. This might not help but just I case ...
> >
> > If you want a closed boundary you need to apply Robin boundary
> conditions because by definition the flux contains two components: one
> related to the diffusion process and one the advection process.
> >
> > For example, http://scicomp.stackexchange.com/a/10576/3691
> >
> > Dan
> >
> > On 24 Aug 2016, at 22:56, Yun Tao <[email protected]> wrote:
> >
> >> Hi all,
> >>
> >> I'm experiencing a bizarre issue when trying to implement zero-flux
> external boundary condition when solving for transient solutions. My
> understanding is that it is the default setting in FiPy 3. Indeed, without
> specifying it, the solutions (attached) remains at unity throughout the
> simulation duration. However, when I manually tried to fix the exterior
> face gradient to zero with var.faceGrad.constrain(0,
> where=m.exteriorFaces), the solution began to blow up (as shown in the
> printed statements). Why does this happen? Is there a hidden conflict in
> conservation settings I should be careful of?
> >>
> >> Thanks,
> >> Yun
> >>
> >>
> >> --
> >> Yun Tao
> >> Postdoc
> >> Center for Infectious Disease Dynamics
> >> Pennsylvania State University
> >> State College, PA 16803
> >> <0flux_bc_fipylistserve.py>
> >> _______________________________________________
> >> 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 ]
> >
> >
> >
> >
> > --
> > Yun Tao
> > Postdoc
> > Center for Infectious Disease Dynamics
> > Pennsylvania State University
> > State College, PA 16803
> > _______________________________________________
> > 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 ]
>
--
Yun Tao
Postdoc
Center for Infectious Disease Dynamics
Pennsylvania State University
State College, PA 16803
from fipy import *
import numpy as np
from fipy import numerix
import matplotlib.pyplot as plt
from matplotlib.mlab import bivariate_normal
import os, sys
fig = plt.figure(figsize=(12, 5))
ax1 = plt.subplot((111))
''' mesh config '''
l = 8.
nx = 101
dx = l/nx
dy=dx
ny = nx
m = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy) + [[-l/2]]
''' set IC; define variables '''
ic = bivariate_normal(m.x, m.y, .1, .1, -3., 0.)
var = CellVariable(mesh=m, value=ic)
''' define system variables & parameters'''
pos = FaceVariable(mesh=m, value=m.getFaceCenters(), rank=1)
xf, yf = pos
D = 1.
convection = VanLeerConvectionTerm
c = 3.
b = ((c,), (c,))
alpha = 2./(dx * .5)
den = (3., 0.)
r = numerix.sqrt((xf-den[0])**2 + (yf-den[1])**2)
faceVelocity = b*numerix.tanh(alpha*r)*(pos-den)/r
''' specify Robin boundary condition '''
#var.faceGrad.constrain([- (1./D) * faceVelocity * var.faceValue], where=m.exteriorFaces)
''' advection-diffusion equation '''
eq3 = TransientTerm() == DiffusionTerm(coeff=D)+ convection(coeff=faceVelocity)
''' viewer config '''
w = 4.
if __name__ == '__main__':
viewer2 = Matplotlib2DGridContourViewer(vars=var,
limits={'xmin': -w, 'xmax': w, 'ymin': -w/2, 'ymax': w/2},
cmap = plt.cm.jet,
axes = ax1)
viewer2.plot()
viewer2.plot()
plt.show()
''' transient analysis '''
dt = 3*.1 * 1./2 * dx / c
steps = 500
O_index_L = []
for step in range(steps):
print 'step ', step
eq3.solve(var=var,
dt = dt)
if __name__ == '__main__':
viewer2.plot()
print 'cell volume: ', var.getCellVolumeAverage() * m.getCellVolumes().sum()_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
