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
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)
pos = FaceVariable(mesh=m, value=m.getFaceCenters(), rank=1)
xf, yf = pos
#var.faceGrad.constrain(0, where=m.exteriorFaces) #weird: adding this line actually blows up the solution
''' Internal BC '''
# middle bar
mask_cc = ((m.x < -2.) | (m.x > 2.) | (m.y < -0.5) | (m.y > 0.5))
mask_ccf = ((xf < -2.) | (xf > 2.) | (yf < -0.5) | (yf > 0.5))
D = 1. * mask_cc
convection = VanLeerConvectionTerm
c = 1.
b = ((c,), (c,)) * mask_ccf
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
eq3 = TransientTerm() == DiffusionTerm(coeff=D)+ convection(coeff=faceVelocity)
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()
dt = 3*.1 * 1./2 * dx / c
steps = 20*np.int(c)
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
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
