FiPy evaluates different solutions on the same system using python2 or python3

I'm experimenting FiPy with a system of partial differential equations and I get different results running the same script in Python2 or Python3. More precisely, I get the result I expect running the script in Python 2.7, while I get a completely wrong result in Python 3.8.

Does somebody knows the reason? Am I missing something? I know that some of the solvers FiPy uses in Python 2.7 don't support Python 3.x, but is this enough to explain a different behaviour?

Thank you in advance.

NOTE: I wanted to post images but I can't since I just subscribed to stack overflow.

Code

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...")

Solution

tl;dr: You're probably seeing a difference in the behavior of the Matplotlib viewer and not in the solver.

When I adjust for bugs and regressions with Matplotlib, I see very similar evolution with Py27 and Py3k. Until we get that fixed, you can see if these changes help:

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)

The solver is different (I get PySparse LinearLUSolver with Py27 and PETSc LinearGMRESSolver with Py3k). This could matter, particularly since you don't sweep this nonlinear problem, but I don't see a huge difference.

The initial condition with Py27 is probably not what you want because of integer division. I recommend initializing with:

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)