Hi FiPy community,
I'm solving a system of 2-dimensional PDEs in time, trying to plot the
changing contour lines of individual solutions together with a map of their
sum values.
In the attached example, I defined four uncoupled diffusion equations with
separate centers. I solved them simultaneously, plotted each solution using
the contour function in matplotlib, and plotted their sum using a FiPy
viewer function.
However, I noticed that the plots are layered only when my FiPy viewer
calls *Matplotlib2DGridContourViewer *(line 55), but not when I try to use
*Matplotlib2DGridViewer *or *Matplotlib2DViewer*, both of which failed to
render any image when the other contour functions are called.
My problem requires me to generate a heat map of the solutions' sum, so I'm
wondering if there's any way to work around this pesky issue.
Thanks,
Yun
--
Yun Tao
Postdoc
Center for Infectious Disease Dynamics
Pennsylvania State University
State College, PA 16803
#!/usr/bin/env python
# encoding: utf-8
"""
Created by Yun Tao on 2012-04-01.
Copyright (c) 2012 __MyCompanyName__. All rights reserved.
"""
from fipy import *
from fipy import numerix
from matplotlib.mlab import bivariate_normal
import matplotlib.pyplot as plt
import numpy as np
fig = plt.figure()
txt = fig.suptitle('time zero', fontsize=14)
plt.show()
ax1 = plt.subplot()
# mesh config
l = 6.
nx = 201.
ny = nx
dx = l / nx
dy = dx
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy) + [[-l/2.]]
# initial condition config
x, y = mesh.getCellCenters()
X = x.reshape(201, -1)
Y = y.reshape(201, -1)
z1 = bivariate_normal(x, y, .2, .2, .3, .3)
z2 = bivariate_normal(x, y, .2, .2, -.3, -.3)
z3 = bivariate_normal(x, y, .2, .2, -.3, .3)
z4 = bivariate_normal(x, y, .2, .2, .3, -.3)
# variable config
phi1 = CellVariable(name='$\phi1$', mesh=mesh, value=z1)
phi2 = CellVariable(name='$\phi2$', mesh=mesh, value=z2)
phi3 = CellVariable(name='$\phi3$', mesh=mesh, value=z3)
phi4 = CellVariable(name='$\phi4$', mesh=mesh, value=z4)
# Boundary Conditions
facesTopRight = ((mesh.getFacesRight())
| (mesh.getFacesTop()))
facesBottomLeft = ((mesh.getFacesLeft())
| (mesh.getFacesBottom()))
BCs = (FixedFlux(faces=mesh.getFacesRight(), value=0.),
FixedFlux(faces=mesh.getFacesLeft(), value=0.))
# Viewer config
if __name__ == '__main__':
viewer = Matplotlib2DGridContourViewer(vars=phi1+phi2+phi3+phi4,
limits={'xmin': -l/4, 'xmax': l/4, 'ymin': -l/4, 'ymax': l/4},
datamin=0.,
axes=ax1,
legend=None)
# equation config
D = 0.1
# define eqs
eq1 = (TransientTerm(var=phi1) == DiffusionTerm(coeff = D, var=phi1))
eq2 = (TransientTerm(var=phi2) == DiffusionTerm(coeff = D, var=phi2))
eq3 = (TransientTerm(var=phi3) == DiffusionTerm(coeff = D, var=phi3))
eq4 = (TransientTerm(var=phi4) == DiffusionTerm(coeff = D, var=phi4))
eq = (eq1 & eq2 & eq3 & eq4)
# temporal config
dt = 0.05
steps = 15
txt.set_text(r'Time Step 0')
for step in range(steps):
print 'steps:', step+1
eq.solve(boundaryConditions=BCs, dt=dt)
txt.set_text(r'Time Step %s'%(step+1))
ax1.cla()
levels=np.arange(0.,.8,.1)
cs1=ax1.contour(X, Y, phi1.value.reshape(201,-1), levels, colors='k', alpha=.3)
cs2=ax1.contour(X, Y, phi2.value.reshape(201,-1), levels, colors='k', alpha=.3)
cs3=ax1.contour(X, Y, phi3.value.reshape(201,-1), levels, colors='k', alpha=.3)
cs4=ax1.contour(X, Y, phi4.value.reshape(201,-1), levels, colors='k', alpha=.3)
ax1.clabel(cs1, fontsize=9, inline=1)
ax1.clabel(cs2, fontsize=9, inline=1)
ax1.clabel(cs3, fontsize=9, inline=1)
ax1.clabel(cs4, fontsize=9, inline=1)
viewer.plot()_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
