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
import numpy as np
from fipy import *
import random
#np.random.seed(13)
from math import *
from scipy import interpolate
import matplotlib.pyplot as plt
import matplotlib.mlab as ml
import time
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 '''
ido = np.random.randint(nxy, size=(100,2)).astype(np.int)
# -- interpolate points
og = np.zeros((nxy, nxy))
np.add.at(og, (ido[:,0], ido[:,1]), 1)
signal_fct = interpolate.interp2d(ulist, ulist, og, kind='cubic')
# signal_fct here is a function that can output the local signal density given any x-y coordinate
''' Convection-Diffusion kernel '''
# -- Parameters
D = 1.
c = 1.
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 = FaceVariable(mesh=m, value=m.getFaceCenters(), rank=1)
xf, yf = m.getFaceCenters()
r = numerix.sqrt((xf-attractor[0])**2 + (yf-attractor[1])**2) # radius of any position to the attractor
faceVelocity = b*numerix.tanh(delta*r)*(s-attractor)/r # convection component in the absence of signal
# -- Fokker-Planck equation
eq = (TransientTerm() == DiffusionTerm(coeff=D)) #+ convection(coeff=faceVelocity * signal.faceValue))
# -- Visualization
if __name__ == '__main__':
viewer = Matplotlib2DGridContourViewer(vars=var,
limits = {'xmin': 0, 'xmax': L, 'ymin': 0, 'ymax': L},
cmap = plt.cm.Greens,
axes = ax1)
''' 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()
ax2.cla()
ax2.plot(ulist[ido[:,1]], ulist[ido[:,0]], 'ko', alpha=.3)
plt.draw()
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
