首页 > 解决方案 > FiPy:我可以根据相邻单元格直接更改 faceVariables 吗?

问题描述

我正在使用二维网格上微生物量 (b1) 分布的生物模型。从生物质中产生蛋白质(p1)。生物质在网格上扩散,而蛋白质则没有。只有产生一定量的蛋白质(p > p_lim),生物量才应该扩散。我尝试通过使用虚拟单元变量 z 乘以扩散系数并仅在 p > p_lim 的单元中将其从 0 设置为 1 来实现这一点。

该条件运行良好,当一个细胞中达到 p 的临界量时,z 设置为 1,并且发生扩散。但是,扩散仍然不能以我想要的速率工作,因为要计算扩散,使用的是面变量,而不是单元格本身的值。z 的面始终是 z=1 的单元格及其 z=0 的相邻单元格的平均值。然而,II 希望扩散以原始速率工作,即使相邻单元仍处于 p < p_lim。

所以,我的问题是:我可以以某种方式访问​​ faceVariable 并更改它吗?例如,如果任何相邻单元已达到 p1 > p_lim,则将面设置为 1?我想这不是一个合适的数学方法,但我想不出另一种方法来模拟这个问题。

我将在下面展示我的模型的一个非常简化的形式。无论如何,我非常感谢您抽出宝贵的时间!

##### produce mesh

nx= 5.
ny= nx
dx = 1.
dy = dx
L = nx*dx
mesh = Grid2D(nx=nx,ny=ny,dx=dx,dy=dy)


#parameters
h1 = 0.5 # production rate of p
Db = 10. # diffusion coeff of b
p_lim=0.1 


# cell variables
z = CellVariable(name="z",mesh=mesh,value=0.)

b1 = CellVariable(name="b1",mesh=mesh,hasOld=True,value=0.)

p1= CellVariable(name="p1",mesh=mesh,hasOld=True,value=0.)


# equations
eqb1 = (TransientTerm(var=b1)== DiffusionTerm(var=b1,coeff=Db*z.arithmeticFaceValue)-ImplicitSourceTerm(var=b1,coeff=h1))
eqp1 = (TransientTerm(var=p1)==ImplicitSourceTerm(var=b1,coeff=h1)) 

# set b1 to 10. in the center of the grid
b1.setValue(10.,where=((x>2.)&(x<3.)&(y>2.)&(y<3.)))
         
vi=Viewer(vars=(b1,p1),FIPY_VIEWER="matplotlib")


eq = eqb1 & eqp1

from builtins import range
for t in range(10):
    b1.updateOld()
    p1.updateOld()
    z.setValue(z + 0.1,where=((p1>=p_lim) & (z < 1.)))
    
    eq.solve(dt=0.1)
     
    vi.plot()

标签: pdefipy

解决方案


除了.arithmeticFaceValue,FiPy 还提供了单元格和面值之间的其他插值器,例如.harmonicFaceValue.minmodFaceValue

这些属性是使用 的子类实现的_CellToFaceVariable,特别_ArithmeticCellToFaceVariable_HarmonicCellToFaceVariable_MinmodCellToFaceVariable

您还可以通过子类化来制作自定义插值器_CellToFaceVariable。两个这样的例子是_LevelSetDiffusionVariableScharfetterGummelFaceVariable(恐怕都没有很好的记录)。

您需要覆盖该_calc_()方法以提供您的自定义计算。此方法接受三个参数:

  • alpha:从面到一侧的单元格的距离与从另一侧的单元格到第一侧的单元格的距离的比率(0-1)的数组
  • id1:面部一侧的单元格索引数组
  • id2:面部另一侧的单元格索引数组

注意:您可以忽略任何子句if inline.doInline:并查看_calc_()该子句下定义的方法else:


推荐阅读