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
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]
''' 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], [[50,0],[0,50]], 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)
''' 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
for step in range(100):
print 'time step', step
eq.solve(var=var, dt = dt)
print 'cell volume: ', var.getCellVolumeAverage() * m.getCellVolumes().sum()
viewer.plot()
viewer2.plot()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]
''' 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], [[50,0],[0,50]], 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], [[50,0],[0,50]], 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
for step in range(100):
print 'time step', step
eq.solve(var=var, dt = dt)
print 'cell volume: ', var.getCellVolumeAverage() * m.getCellVolumes().sum()
varS = update_signal(step/3.)
print 'max signal location', m.getCellCenters()[:, numerix.argmax(varS)] # the overall increase in the first element shows that signals are indeed getting updated
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 ]
