How can one apply "vectorized=True" to the Scipy solver?

[HannoSpreeuw]

Take, e.g. the SIR example.
Replace this line

sol = eq.solve(state, t_range=50, dt=1e-2, tracker=["progress", tracker])

by these

solver = ScipySolver(eq)
controller = Controller(solver, t_range=50, tracker=["progress", tracker])
sol = controller.run(state)

and it runs just as smooth.
However, when I change that first line to

solver = ScipySolver(eq, method = "Radau", vectorized = True)

I get

ValueError: cannot reshape array of size 9437184 into shape (3,32,32)

I am puzzled, I am not sure in what way I should change evolution_rate to make this work. Or do I need to change something else?

[david-zwicker]

I don't think vectorize=True will ever be supported. Have you tried without this option?

[HannoSpreeuw]

Perhaps I misinterpreted this comment

This class is a thin wrapper around :func:`scipy.integrate.solve_ivp`. In
    particular, it supports all the methods implemented by this function.

from the pde.solvers.scipy source code.

On the other hand, it should not be too hard to make Jacobians work within py-pde, right?
Just requires an extra helper for the Jacobian, similar to rhs_helper for fun, the right-hand side of the system.

[HannoSpreeuw]

But thanks for the tip about tweaking the tolerances for the explicit solver.

For my set of pdes I used

solver = ExplicitSolver(eq, scheme="rk", adaptive=True, backend="numba", tolerance=1e-3) 

with tolerances between 1e-2 and 1e-6 and this always resulted in

RuntimeError: Time step below 1e-10

So these equations must be stiff.

[gnorman7]

I'm confused by this thread. The source code for the ScipySolver has the following, within the stepper method:

sol = integrate.solve_ivp(
    rhs_helper,
    t_span=(t_start, t_end),
    y0=np.ravel(state.data),
    t_eval=[t_end],  # only store necessary data of the final time point
    **self.solver_params,
)

And further, calling solver.solver_params does show that the kwargs are being passed properly. Further, trying the input method='foo' gives a SciPy error that 'foo' is not recognized as a solver name. So it seems like it is set up to integrate with SciPy. However, I'm not able to get different behavior by switching the method or the tolerances.

[david-zwicker]

What do you mean by "I'm not able to get different behavior by switching the method or the tolerances.". What have you tried and what did you expect?

[gnorman7]

Nevermind that part, I realize now that I had an error in my code.