Issue with correctness of the results

[ma-sadeghi]

I was trying to integrate the heat equation on a 100x100x100 grid, and scikits.odes's rk5 gives incorrect results (compared to scipy's solve_ivp). The values I get from scipy are always bounded between 0 and 1 (which they should be since I use 0 and 1 as Dirichlet boundary conditions), but those I get from scikits.odes blow up to ~10^10. Here's how to reproduce it:

# %% Import libraries
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
from scikits.odes.odeint import odeint
from time import time

# %% Define du/dt
def func(t, v, du, params):
    h, k, Nx, Ny, Nz = params
    is3d = True if Nz > 1 else False
    u = v.reshape((Nx, Ny, Nz))
    idx = np.arange(1, Nx-1)
    idy = np.arange(1, Ny-1)
    idz = np.arange(1, Nz-1) if is3d else np.array([0])
    u_xx = (u[np.ix_(idx-1, idy, idz)] - 2*u[np.ix_(idx, idy, idz)] + u[np.ix_(idx+1, idy, idz)]) / h**2
    u_yy = (u[np.ix_(idx, idy-1, idz)] - 2*u[np.ix_(idx, idy, idz)] + u[np.ix_(idx, idy+1, idz)]) / h**2
    if is3d:
        u_zz = (u[np.ix_(idx, idy, idz-1)] - 2*u[np.ix_(idx, idy, idz)] + u[np.ix_(idx, idy, idz+1)]) / h**2
    else:
        u_zz = 0.0
    idx1d = np.ravel_multi_index(np.ix_(idx, idy, idz), (Nx, Ny, Nz)).flatten()
    du[idx1d] = k * (u_xx + u_yy + u_zz).flatten()
    return du

# %% Define model parameters
Nx, Ny, Nz = 100, 100, 100
Lx = 1.0
Ly = Lx * (Ny-1)/(Nx-1)
Lz = Lx * (Nz-1)/(Nx-1)
h = Lx/(Nx-1)
D = 5e-3
params = (h, D, Nx, Ny, Nz)

# Safety checks
assert Nx >= 3
assert Ny >= 3
if Nz != 1:
    assert Nz >= 3

# %% Scikits.Odes
def func_odes(t, u, du):
    func(t, u, du, params)

tspan = (0.0, 5.0)
tout = np.linspace(tspan[0], tspan[1], 10)
u0 = np.zeros((Nx, Ny, Nz))     # IC @ t = 0 ; u = 0.0
u0[0, ...] = 1.0                # BC @ x = 0 ; u = 1.0
u0[-1, ...] = 1.0               # BC @ x = L ; u = 1.0
u0[:, 0, :] = 1.0               # BC @ y = 0 ; u = 1.0
u0[:, -1, :] = 1.0              # BC @ y = L ; u = 1.0

t0 = time()
out = odeint(func_odes, tout, u0.flatten(), method="rk5",
             rtol=1e-3, atol=1e-3)
print(f"Scikits.Odes runtime (rk5): {time()-t0:.2f} s")

# t0 = time()
# out = odeint(func_odes, tout, u0.flatten(), method="cvode", rtol=1e-3, atol=1e-3)
# print(f"Scikits.Odes runtime (cvode): {time()-t0:.2f} s")

u3d_odes = out.values.y.reshape((tout.size, Nx, Ny, Nz))

fig, ax = plt.subplots()
f = ax.imshow(u3d_odes[-1, ..., Nz//2], vmin=0, vmax=u3d_odes.max())
ax.set_title("Scikits.Odes")
plt.colorbar(f)

# %% Scipy
du = np.zeros(u0.size)
def func_scipy(t, u):
    return func(t, u, du, params)

t0 = time()
sol = solve_ivp(func_scipy, tspan, u0.flatten(), method="RK45",
                t_eval=tout, rtol=1e-3, atol=1e-3)
print(f"Scipy runtime (rk45): {time()-t0:.2f} s")

u3d_scipy = sol.y.reshape(Nx, Ny, Nz, sol.t.size)

fig, ax = plt.subplots()
f = ax.imshow(u3d_scipy[..., Nz//2, -1], vmin=0, vmax=u3d_scipy.max())
ax.set_title("Scipy")
plt.colorbar(f)
plt.show()

[bmcage]

This is not a problem with odes, but with the solver rk45, which is in the background the dopri5, which is actually from scipy.integrate also, see
https://github.com/bmcage/odes/blob/master/scikits/odes/dopri5.py#L104

So you should be able to reproduce this with scipy, and raise a ticket there. However, as this will be the underlying solver, not something scipy will actually solve, the solution is always to use the solver as intended, which you are not doing. Looking at the scipy code, odeint rk5 and and RK45 of solve_ivp are probably the same...

A good way forward is to use the step methods, and keep track of steps, and errors, along the way.

[bmcage]

You made the error of assuming in odes that du is initialized 0, or copy of before, which it is not. You could use own memory for this, but for now, simplest solution is adding init to 0:

du[:] = 0
du[idx1d] = k * (u_xx + u_yy + u_zz).flatten()
return du

Like this rk5 of odes will give you the correct result.
You could deduce this from noting that the only reason the border can change value is if du is not 0 on the border. With dirichlet conditions, not including them in the system is a saver bet.

[ma-sadeghi]

Again, thank you very much for your thorough response :)