Add more (default) methods to problem class interface

[tlunet]

Following a post by @brownbaerchen in #244 (comment) , the current interface for Problem class in pySDC is as follow :

class ptype(object):
    def __init__(self, init, dtype_u, dtype_f):
        # ...
     def eval_f(self, u, t):
        # ABSTRACT
    def apply_mass_matrix(self, u):
        # ABSTRACT

Looking at different current problem implementations, here are some ideas

• define by default apply_mass_matrix(self, u) as an identity operator on u (since most problem don't have a mass matrix) and suggest to redefine it in a children class using, for instance, FEM.
• add u_exact(self, t, u_init=None, t_init=None) as a standard method using the scipy reference solution for IVP, and allow it to be redefined for problem that have a simpler exact solution
• add solve_system(self, rhs, factor, u0, t) as a base method to the problem class, with a generic implementation of a solution of the (non)-linear system u-factor*f(u,t) = rhs using the (very generic) Newton-Krylov solver of scipy, allow it to be redefined for problems that have a simpler way to solve this (e.g for linear problems) or even for problem that don't allow any implicit solve (and use a custom ImplicitNotAllowedError).

📜 The change on apply_mass_matrix and solve_system would facilitate the binding with dedalus V3, as those are two main components required by dedalus time-integrators ...

Any comment/suggestion on that @pancetta @brownbaerchen ?

[pancetta]

When we first added FEniCS we needed to rewrite the sweeper and transfer classes, because now we had a mass matrix. If you don't have a mass matrix, those classes are easier to read and digest (and much closer to the literature). This is why we decided against a default apply_mass_matrix function. So far, only FEM needs that and this is what we're doing only 1% of the time. But maybe we should actually do that more often anyway?

[tlunet]

So far from what I saw in dedalus, they also have many situations where the mass matrix is actually identity. But their implementation solves generic PDE problem of the form (or something like)

$$M \frac{dU}{dt} = L U + f(U,t)$$

with $L$ for the linear term and f(U,t) for the non linear term (solved with IMEX). Considering $M=I$ and $L=0$ allows to retrieve the classical form of ODE systems. Maybe a similar approach could be considered in pySDC, such that the problem interface itself allows complex ODE system formulations, but the default implementation sticks to the classical ODE representation ? That is :

$$\frac{dU}{dt} = f(U,t)$$

[pancetta]

Sure, makes sense. One additional case where a mass matrix is needed is for compact FD stencils.. so, yes, maybe it's about time we don't treat this as an edge case but rather as a more generic way to work with PDE discretizations.

[tlunet]

After thinking about it on a (long ...) train travel, here are some propositions.

Mass-matrix

Add it by default to define the problem

$$M\frac{dU}{dt} = f(U,t)$$

and define a base method applyMassMatrix(u) that would simply return u in the base class. If someone wants a non-identity mass-matrix (e.g for FEM or compact FD), its should simply be redefined in the specialized problem class. It also means that the generic sweeper should later use this applyMassMatrix method by default.

IMEX

Don't split $f(U,t)$ into non-linear and linear terms, for instance

$$M \frac{dU}{dt} = L U + F(U,t)$$

since some application may consider applying implicit sweep to a non-linear term, or a time-dependent $L(t)$ matrix.
But keep the current generic ODE form with mass matrix (i.e generic DAE form) and use the current idea of the data_f type, that can contains an impl and expl attribute to separate both IMEX part. Maybe a base class should be defined for pySDC datatypes in core module, with the IMEX datatype class (inheriting from base datatype) also defined in this core module.

Additionally, the IMEX sweep could be added to the core modules, close to the base sweeper class that would contains both explicit and implicit implementations (and just testing the QDelta matrix to choose which one to use). That way, for one given application, only datatype and problem classes implementation is required, and the rest is only the correct choice of core classes and parameters already provided by pySDC core.

Generic solver for implicit sweep

Add a generic solveSystem(rhs, factor, u0, t) that use the Krylov-Newton solver from Scipy to solve

$$U - factor*f(U,t) = rhs$$

with u0 as starting solution for the iterative solver. This method could be redefined for linear problems or when the Jacobian matrix of f(U,t) is available and allow a more efficient solve (e.g for the Lorenz system).

Generic ODE solver for u_exact

Define a base generic method uExact(t, u0=None, t0=0) which would use the Scipy approach from @brownbaerchen to solve the ODE system when the applyMassMatrix() method is not redefined (quick easy test to do in Python). If u0 is None, uses the initial solution associated to the problem, given by a u0 property set as abstract in the base class, that must be given with each problem implementation.

For cases when the applyMassMatrix function is redefined, use a more complex approach for generic DAE (see for instance this stackoverflow post).

Would you validate/correct that @pancetta ?

[pancetta]

Thanks a lot @tlunet! Let me try to comment one by one:

• Mass-matrix
  Yes! We have to remember that this changes the base transfer, too.

• IMEX
  pySDC does not enforce linear/nonlinear splitting. It's up to the user to define what f.expl, f.impl or solve_system means. I'd rather not move the sweeper implementation to the core classes, this is a clear implementation thing. We can think about setting defaults, yes, but from the implementation classes.

• Generic solver for implicit sweep
  Interesting thought, I like it!

• Generic ODE solver for u_exact
  First part: yes. Second part: We're currently adding DAE support to pySDC (through a Master's project) and I don't want to interfere here. Also, I think only the most simplest things should be set by default and if users diverge, they need to take care of the implications themselves. If it turns out to be necessary (depending how people make use of the mass matrix feature), we can still reconsider.

[brownbaerchen]

I have no idea about mass matrices, so I can't give any comments on that.

I do however like the idea of implementing IMEX on a more basic level. In particular, I would like to have a function like

eval_f_fully_implicit

which would in the case of IMEX return the sum of the two rhs evaluations and in the implicit sweeper just the implicit part and so on.
Since I never used any other sweeper, I don't know if this can be generalised to those as well.

This would allow to have a bit more abstraction in the sweeper classes. For instance, you wouldn't need to overload the integrate function in the IMEX sweeper anymore.

The other place where this is useful is the scipy stuff. There you need to supply something that evaluates the full rhs. Ofc naming it "fully_implicit" may not make sense, since in this case "fully_explicit" is probably closer to the truth, but whatever.

We can also just implement a function like this for the scipy reference solution stuff, if @pancetta wants to keep the non-imex sweepers a bit more simple, though.

[tlunet]

I do like this idea of full rhs evaluation, but I wonder if this should be done rather in the data_f class than the problem class.

I would imagine that the base data_f class must have a vals attribute that would expose data values in the generic case, e.g the numpy array for mesh datatype, and for IMEX datatype it would simply sum the implicit and explicit part and return it.

Eventually, this could also be a property that would prevent writing on it, except if you are in the problem class, and then fully aware of what the data_f type is, and then using the specialized data setter (e.g f.vals[:,:,:] = ... for a 3D problem)

[pancetta]

I spend a lot of time getting rid of a dedicated attribute (called values in previous versions, pre-v4), because then you need to write functions to do stuff like add, rmult and so on (or override those methos). This usually involves instantiating new data types each time you do that and this is very costly. If you can avoid this, fine, if not, I won't go back to dedicated attributes.

[pancetta]

The way it's currently done/preferred is to subclass an existing datatype (e.g. ndarray) and modify the __new__ function. This way all operations are preserved and no further overhead is added.

[tlunet]

Should I consider then that the fenics_mesh class is supposed to be deprecated ?

[pancetta]

Not all data structures can be subclassed. The one from FEniCS are tricky. I'm trying to get FEniCSx to work, there I'd probably rely on mesh as a wrapper. So, no, not deprecated.

[tlunet]

So ... it appears that implementing those generic methods for problem classes need also a rework of the datatype classes, with the introduction of a base interface for u and a base generic class for evalF output, that would be directly linked to the datatype of u.

Here is a first draft that tries to solve those issues : https://github.com/Parallel-in-Time/pySDC/blob/847e8ca1491df1076aff9c61768245ff940101a9/pySDC/playgrounds/datatypes/draft.py. I'm slowly checking how the current class could be translated and used into this new structure. Let me know if you see any issue and supplementary test I could try.

@pancetta : for some reason I don't understand yet, it seems that we don't need the complicated approach to subclass ndarray, as described here ...

[tlunet]

PPS : we don't want to do arithmetic operations between DataTypeF types, right ?

What do you mean? Adding two different data types (e.g. fenics_mesh and cupy_mesh)? No. We don't want that.

No, I mean adding the result of two eval_f evaluation, for instance. In general, for one given application, adding two variables that have the same dtype_f type ...

[pancetta]

I'm confused and pretty sure I don't understand what you're saying: f1 + f2 is a common operation, e.g. during integration, no? So, yes, of course we need to add f-evaluations.

[tlunet]

I mean : you never dtype_f + dtype_f, but rather dtype_f + dtype_u, that you store into dtype_u, right ?

[pancetta]

Sure you do dtype_f + dtype_f! Computing Qf needs exactly that.

[pancetta]

Oh, I see... hm. Very good question...

[pancetta]

Can you please test parallel code, too? E.g. distributed arrays with a communicator, computing a norm or so.

[pancetta]

So, we need different communicators for different instances of data. One u could belong to this space communicator, another u to another one. Can you do that with your approach?

[pancetta]

Appealing II, I must admit. But again, we need to check for performance.

[tlunet]

Still ongoing ... but here are some preliminary results with one Forward Euler step. Dashed lines are the current implementation (reference) with mesh datatypes.
The u + dt*f(u,t) test consider both evaluation of f and multiplication, while the u + dt*f test only look at the FE update with the f(u,t) term already computed. See the current test file on v6.

[pancetta]

So, you're telling me your approach may even be faster...?

[tlunet]

Not promising anything yet, but seems like it 😃

[pancetta]

We need to test this with PETSc and FEniCS, too. And I worry a bit about the GPU stuff...

[pancetta]

Bottom line: I don't see the benefit of rewriting the data types yet again. If the way you propose it faster, I'm on board. If not, I'd rather leave it like that, because we have so many open issues, even without touching the data types.

[tlunet]

This would allow :

• to write a generic solveSystem and uExact on the base level
• one unique implementation of integrate that would work for any datatypes
• one unique implementation of a generic sweeper (that would call solveSystem only for non-null diagonal QDelta coefficients),
• an easy definition of IMEX datatypes (or multicomponents for multi-implicit SDC) that could be combined with any main datatypes (mesh, petsc, fenics, particles, ...)
• ...

So mostly bonus things that may be not necessary yet, but could surely improve usage and development of pySDC. But I do agree : this must be adapted to the current datatype, and at least not slower than the current implementation (better if faster).

[pancetta]

Thanks! Don't me wrong, I really appreciate your effort and dedication!

to write a generic solveSystem and uExact on the base level

This depends on the actual data type anyway, no? solveSystem is very different for FEniCS than for PETSc.

one unique implementation of integrate that would work for any datatypes

Ok, well, but this is going to less explicit, right? In general the only dependency is on the number and type of f components. The rest (which specific type) is generic.

one unique implementation of a generic sweeper (that would call solveSystem only for non-null diagonal QDelta coefficients)

Which sweepers would it replace? IMEX? implicit? explicit? Verlet? Multi-implicit?

an easy definition of IMEX datatypes (or multicomponents for multi-implicit SDC) that could be combined with any main datatypes (mesh, petsc, fenics, particles, ...)

Appealing, indeed!