Issue in specifying explicit space dependence in 1D PDE custom class

[Raffo921]

Hello to all the community members!

First things first, I gratefully thank you to all the creators and maintainers for your hard work in creating this amazing library. It is an invaluable tool with tons of functionalities for those involved in Scientific computing.

I have been using it for few weeks so I'm quite new to its use. My goal now is to implement a custom class PDE which has explicit time and space dependence. Instead of jumping directly to my problem, I tried to modify some examples in the Examples to see how the solution behaves.

In the Korteweg-de Vries equation, I modified the evolution_rate by simply adding an explicit x dependence as state.grid.axes_coords[0], taking a leaf from Q&A #566.

def evolution_rate(self, state, t=0):
        """Implement the python version of the evolution equation."""
        assert state.grid.dim == 1  # ensure the state is one-dimensional
        
        grad_x = state.gradient("auto_periodic_neumann")[0]
        return 6 * state*state.grid.axes_coords[0]*grad_x - grad_x.laplace("auto_periodic_neumann")

The problem is that it returns this error AttributeError: 'list' object has no attribute 'reshape' . If I use instead state.grid.cell_coords[0] it works, giving this field

(which I don't know if it is the correct solution to this equation)

I tried the same thing with 2D PDEs by modifying the FitzHugh–Nagumo example with expression state.grid.axes_coords[0]

def evolution_rate(self, state, t=0):
        v, w = state  # membrane potential and recovery variable

        v_t = v.laplace(bc=self.bc) + v*state.grid.axes_coords[0] - v**3 / 3 - w + self.stimulus
        w_t = (v + self.a - self.b * w) / self.τ

        return FieldCollection([v_t, w_t])

and in this case it works, giving

Explicit time dependence seems working in both 1D and 2D PDEs. I simply used t for time in the return of evolution_rate since it is an argument.

My question is, Am I missing something or it could be a bug in case of 1D problems?
Thank you in advance

[david-zwicker]

Thank you for your detailed description of the problem, which helped me to re-create it on my side. It turns out that the problem you observed was a combination of an unstable equation and a suppressed error message. You thus used the package correctly, but the equation you describe blows up heavily, which py-pde cannot deal with. The error is somewhat more transparent when you choose solver="explicit", and I just modified the package to also emit a more useful error when the scipy solver fails (see #679).

[Raffo921]

Dear David,

Thank you for your fast reply.

I thought about instability, but I expected a different behavior. Like the solver ran smoothly without errors or overflows, but the solution magnitude would be like 10^30.

[david-zwicker]

Yes, I agree that the error messages could be more helpful. With the latest version of our package, we are now at least exposing the error message of scipy, but this is still basically only saying that the step failed. I guess the issue is that scipy uses adaptive time stepping, so the solution does not explode, but rather it becomes impossible to find a suitable time step.