FiPy 使用 python2 或 python3 代码

问题描述

我正在用偏微分方程组试验 FiPy，我在 Python2 或 python3 中运行相同的脚本得到不同的结果。更准确地说，我得到了我期望在 Python 2.7 中运行脚本的结果，而我在 Python 3.8 中得到了完全错误的结果。

有人知道原因吗？我错过了什么吗？ 我知道 FiPy 在 Python 2.7 中使用的一些求解器不支持 Python 3.x，但这足以解释不同的行为吗？

提前致谢。

注意：我想发布图片，但我不能，因为我刚刚订阅了堆栈溢出。

代码

import fipy as fp
import fipy.tools as fpt

"""
Let's try to solve the reaction diffusion equation fith FiPy
    A + B ---k---> C
    
"""
# define spatial domain (square)
nx = ny = 20
dx = dy = 1.
mesh = fp.Grid2D(dx=dx,dy=dy,nx=nx,ny=ny)
linear_dimension = (nx * ny,)

# define molecule as CellVariables
a = fp.CellVariable(name="molecule A",mesh=mesh)
b = fp.CellVariable(name="molecule B",mesh=mesh)
c = fp.CellVariable(name="molecule C",mesh=mesh)


# define initial condition
def gauss_2d(mu_x,mu_y,sigma_x,sigma_y,x_1d,y_1d):
    """
    Utility function to define initial conditions (see below). Provide a simil-gaussian
    distribution (NOTICE: it's not an exact gaussian because I needed it to be 1 at the top) 
    """
    # initialize matrix
    gauss_mat = fpt.numerix.empty((len(x_1d),len(y_1d)))
    # define gaussian
    gauss_x = fpt.numerix.exp((-1 / 2) * (((x_1d - mu_x) ** 2) / (sigma_x ** 2)))
    gauss_y = fpt.numerix.exp((-1 / 2) * (((y_1d - mu_y) ** 2) / (sigma_y ** 2)))
    # evaluate each point of the matrix
    for i in range(0,len(x_1d)):
        for j in range(0,len(y_1d)):
            gauss_mat[i,j] = gauss_x[i] * gauss_y[j]

    return gauss_mat


normal_distribution = gauss_2d(mu_x=nx / 2,mu_y=ny / 2,sigma_x=fpt.numerix.sqrt(10),sigma_y=fpt.numerix.sqrt(10),x_1d=fpt.numerix.arange(0,nx,dx),y_1d=fpt.numerix.arange(0,ny,dy))
a_max = 100.
a0 = a_max * normal_distribution
a.setValue(a0.reshape(linear_dimension))
b_max = a_max / 2
b0 = b_max * normal_distribution
b.setValue(b0.reshape(linear_dimension))
c0 = 0.
c.setValue(c0)

# create viewer for the three molecules
vi = fp.Viewer(vars=(a,b,c),datamin=0.,datamax=a_max)
vi.plot()
fp.input("Press enter to continue...")

# define the reaction term
k = 0.01
reaction_term = k * a * b

# define the equation for each molecule
D = 0.05
eq_a = fp.TransientTerm(var=a) == fp.DiffusionTerm(coeff=D,var=a) - reaction_term
eq_b = fp.TransientTerm(var=b) == fp.DiffusionTerm(coeff=D,var=b) - reaction_term
eq_c = fp.TransientTerm(var=c) == fp.DiffusionTerm(coeff=D,var=c) + reaction_term

# couple equations
sys = eq_a & eq_b & eq_c

# solve
dt = 0.25
steps = 400
for step in range(steps):
    sys.solve(dt=dt)
    vi.plot()

fp.input("Press enter to continue...")

解决方法

tl;dr：您可能会看到 Matplotlib 查看器的行为与求解器的行为不同。

• 当我使用 Matplotlib 调整错误和回归时，我看到 Py27 和 Py3k 的演变非常相似。在我们解决这个问题之前，您可以看看这些更改是否有帮助：

  diff --git a/rxndiff.py b/rxndiff.py
  index b41932a..ea30a95 100644
  --- a/rxndiff.py
  +++ b/rxndiff.py
  @@ -1,5 +1,6 @@
  import fipy as fp
  import fipy.tools as fpt
  +from matplotlib import pyplot as plt
  
  """
  Let's try to solve the reaction diffusion equation fith FiPy
  @@ -79,5 +80,8 @@ steps = 400
  for step in range(steps):
       sys.solve(dt=dt)
       vi.plot()
  +    for vw in vi.viewers:
  +        vw.fig.canvas.draw_idle()
  +    plt.show(block=False)
  

• 求解器不同（我用 Py27 得到 PySparse LinearLUSolver，用 Py3k 得到 PETSc LinearGMRESSolver）。这可能很重要，特别是因为您没有解决这个非线性问题，但我看不出有什么大的不同。
• Py27 的初始条件可能不是您想要的，因为整数除法。我建议初始化：

  mu_x=nx / 2.
  mu_y=ny / 2.
  sigma_x=fpt.numerix.sqrt(10.)
  sigma_y=fpt.numerix.sqrt(10.)
  
  normal_distribution = (fpt.numerix.exp((-1./2) * (mesh.x - mu_x) ** 2 / sigma_x ** 2)
                        * fpt.numerix.exp((-1./2) * (mesh.y - mu_y) ** 2 / sigma_y ** 2))
  a_max = 100.
  a0 = a_max * normal_distribution
  a.setValue(a0)
  b_max = a_max / 2
  b0 = b_max * normal_distribution
  b.setValue(b0)
  c0 = 0.
  c.setValue(c0)