Hi Jonathan,
I finally got a chance to test out the solution. However, it still doesn't
work and I suspect I'm overlooking something relatively. I have made the
following changes:
- inserted *time*=Variable(): line 25
- updated the time-dependent variable *varS* using *time* in the while
loop: line 101
- updated time.value during the iterations: line 103
I still do not see changes in either *var *or *varS*. My modified code is
attached below -- note that the time-dependency doesn't operate according
to a source term but rather relates to the component inside the convection
term: *varS.faceValue*.
Is my use of an embedded function *update_signal* inside the iteration loop
messing things up?
Thanks,
Yun
On Mon, Oct 19, 2015 at 4:35 PM, Guyer, Jonathan E. Dr. <
[email protected]> wrote:
> You would do this in a similar way to how time-dependent boundary
> conditions are illustrated in examples/diffusion/mesh1D.py:
>
> >>> phi = CellVariable(mesh=..., value=...)
> >>> time = Variable()
> >>> source1 = phi * time
>
> or, e.g.,
>
> >>> source2 = mesh.x * time
>
>
> Then as you iterate in timesteps, you would update time with:
>
> >>> time.value = time.value + dt
>
> source1 and source2 will automatically reflect the new value of time.
>
> On Oct 12, 2015, at 4:15 PM, Yun Tao <[email protected]> wrote:
>
> > Hi FiPy community,
> >
> > This is a potentially helpful update to a recent question I submitted,
> where my goal was to spatially vary the convection strength of a
> Fokker-Planck equation for random variable var(x,t) as a simple function of
> local signal distribution varS(x,t). Specifically, I wanted to solve for
> the transient dynamics of var(x,t), given that, at each location x, its
> probability surface is pulled towards a fixed, central point-attractor to a
> degree that is proportional to the estimated value of varS(x,t).
> >
> > I like to thank Jon for his tremendous amount of help, from which I was
> able to generate a working script on the condition that varS(x,t) is
> time-invariant. The code is attached here as signals_static.py. var(x,t) is
> plotted over time in the left panel of the animation, and varS(x) is on the
> right. The increasing topological distortion shown in the simulation is
> consistent with how we expect var(x,t) to behave.
> >
> > My current goal is to solve for var(x,t) on the condition that varS(x,t)
> also varies, partially stochastically, over time. Note that this doesn't
> involve coupling the Fokker-Plank equation with another differential
> equation. I've attempted to do this in the attached script:
> signals_dynamic.py. The only addition from the previous script is a "signal
> updating function', through which we force the spatial distribution of
> varS(x,t) to shift rightward every time step. However, for some reason,
> var(x,t) is unresponsive to these changes.
> >
> > Therefore, my question is: how can I base the convection term on a
> CellVariable that gets temporally updated outside of the equation
> definition?
> >
> > Thanks,
> > Yun
> >
> > On Wed, Aug 26, 2015 at 2:47 PM, Guyer, Jonathan E. Dr. <
> [email protected]> wrote:
> > Yun -
> >
> > I've gotten your script to "work" and posted the changes to:
> >
> > https://gist.github.com/guyer/caca956463dfc3835722/revisions
> >
> > The main changes I made were:
> >
> > * to get rid of the intrep2d, as it wasn't working properly
> [signal_fct(xf, yf) generates a result of shape
> > (len(xf), len(xf)) instead of (len(xf),).] I was able to get it
> working a bit better, but not completely, and I
> > realized that it doesn't really do anything for you that simply
> placing your signals in a CellVariable and then
> > letting it calculate its .faceValue doesn't accomplish.
> >
> > * simplify the calculation of faceVelocity (m.faceValues is already a
> rank-1 FaceVariable)
> >
> > Although this script functions, I suspect it's not really what you're
> looking for. The signals are all extremely localized and faceVelocity is
> really not responsive to the density of signals, but just discretely to
> whether there's a signal in a given cell. If that's so, I think you'll want
> to calculate a density field for the signals, rather than placing them in
> discrete locations.
> >
> > - Jon
> >
> > On Aug 19, 2015, at 8:00 PM, Yun Tao <[email protected]> wrote:
> >
> > > Hi FiPy community,
> > >
> > > I'm currently trying to combine the powerful tool of FiPy with
> agent-based modeling. The problem I'm trying to solve is this:
> > >
> > > In a 2D landscape scattered with "deterrent point signals", I want to
> solve for the transient solution of a convection-diffusion (Fokker-Planck)
> equation that increases its advection towards its central attractor in a
> way that is proportional to the interpolated density of local signals. I
> therefore expect to see gradual deformation, and slowing down of spread, in
> the solution boundary as diffusion brings it closer to clustered signals.
> > >
> > > However, since the point signals are located on mesh cell centers and
> the convection coefficient in FiPy requires FaceVariable inputs, there is a
> problem with dimensionality I cannot quite understand. How should I
> integrate these two processes?
> > >
> > > I've attached my current script, which has the convection term
> commented out for now. Left figure is the PDE solution; right figure is the
> locations of the signal points.
> > >
> > > Any help would be greatly appreciated.
> > >
> > > Thanks,
> > > Yun
> > >
> > >
> > > --
> > > Yun Tao
> > > PhD
> > > University of California, Davis
> > > Department of Environmental Science and Policy
> > > One Shields Avenue
> > > Davis, CA 95616
> > >
> <fipy-ibm2_test_forum.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
> > PhD
> > University of California, Davis
> > Department of Environmental Science and Policy
> > One Shields Avenue
> > Davis, CA 95616
> >
> <signals_static.py><signals_dynamic.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
PhD
University of California, Davis
Department of Environmental Science and Policy
One Shields Avenue
Davis, CA 95616
import numpy as np
from fipy import *
import random
np.random.seed(13)
from math import *
from scipy.signal import fftconvolve
from scipy import interpolate
import matplotlib.pyplot as plt
import matplotlib.mlab as ml
fig = plt.figure(figsize=(15,6))
ax1 = plt.subplot((121))
ax2 = plt.subplot((122))
plt.ion()
plt.show()
''' Landscape parameters '''
L = 10.
nxy = 100.
dxy = L/nxy
m = Grid2D(nx=nxy, ny=nxy, dx=dxy, dy=dxy)
ulist = m.x.value[:nxy]
time = Variable()
''' Deterrent signals'''
# -- individual signal kernel (radial basis)
bandwidth = .05
i, j = ulist.reshape(nxy,1), ulist.reshape(1,nxy)
r = np.minimum(i-ulist[0], L-i+ulist[0])**2 + np.minimum(j-ulist[0], L-j+ulist[0])**2
rbf = sqrt(1 / (2 * bandwidth ** 2))
ker = np.exp(-(rbf * r) ** 2)
ker = ker/np.sum(ker) **.5
# -- moving the kernel to the center of the lattice in order to apply sp.fftconvolve
ker = np.roll(ker, np.int(nxy//2-1), 0)
ker = np.roll(ker, np.int(nxy//2-1), 1)
# -- randomize signals' distribution
ido = np.random.multivariate_normal([55,35], [[0,0],[0,0]], 100)
ido = np.floor(ido).astype(int)
# -- convolve & estimate the accumulative signal kernel
og = np.zeros((nxy, nxy))
np.add.at(og, (ido[:,0], ido[:,1]), 1)
dd = fftconvolve(og, ker, mode='same')
varS = CellVariable(mesh=m, value=dd.flat)
''' Signal update function '''
def update_signal(rate):
ido = np.random.multivariate_normal([55,35+rate], [[0,0],[0,0]], 100)
ido = np.floor(ido).astype(int)
og = np.zeros((nxy, nxy))
np.add.at(og, (ido[:,0], ido[:,1]), 1)
dd = fftconvolve(og, ker, mode='same')
varS = CellVariable(mesh=m, value=dd.flat)
return varS
''' Convection-Diffusion kernel of `var` '''
# -- Parameters
D = 3.
c = 3.
b = ((c,),(c,))
delta = 2./(dxy * 0.5) # smoothing parameter
convection = VanLeerConvectionTerm
# -- Initialization
attractor = (3., 5.)
z = ml.bivariate_normal(m.x, m.y, .1, .1, attractor[0], attractor[1])
var = CellVariable(mesh=m, value=z)
# -- Spatial objects assignments
s = m.faceCenters
r = (s-attractor).mag
faceVelocity = c*numerix.tanh(delta*r)*(s-attractor)/r
# -- Fokker-Planck equation
eq = (TransientTerm() == DiffusionTerm(coeff=D) + convection(coeff=faceVelocity * varS.faceValue))
# -- Visualization
if __name__ == '__main__':
viewer = Matplotlib2DGridContourViewer(vars=var,
limits = {'xmin': 0, 'xmax': L, 'ymin': 0, 'ymax': L},
cmap = plt.cm.Greens,
axes = ax1)
viewer2 = Matplotlib2DGridContourViewer(vars=varS,
limits = {'xmin': 0, 'xmax': L, 'ymin': 0, 'ymax': L},
cmap = plt.cm.Reds,
axes = ax2)
''' Simulation '''
dt = 0.1 * 1./2 * dxy
steps = 30
while time() < steps * dt:
print 'time step', time
eq.solve(var=var, dt = dt)
print 'cell volume: ', var.getCellVolumeAverage() * m.getCellVolumes().sum()
varS = update_signal(time*50)
print 'max signal location', m.getCellCenters()[:, numerix.argmax(varS)] # the overall increase in the first element shows that signals are indeed getting updated
time.value = time.value + dt
# print var.value[1000] # evaluate value at random position to see if the result changes
viewer.plot()
viewer2.plot()
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
