Implementing other DAE solvers in `pySDC`

[lisawim]

Since the current implemented fully-implicit SDC sweeper for DAEs fully_implicit_DAE and its semi-explicit sibling SemiExplicitDAE should also be competitive against other DAE solvers, sweepers such as Forward Euler, Backward Euler, Runge-Kutta and Multistep shall be implemented. A template is already provided by the Runge-Kutta and Multistep sweepers.

First of all, their framework has to be checked how they do work. @brownbaerchen already adapted the classes in the way, that only the update_nodes method has to be overwritten in case of inheritance from these ones. Since the pull request #371 resulted in a deeper look into the module pySDC.core.Sweeper and how does an inheritance of fully_implicit_DAE from generic_implicit could work, also the methods compute_residual and predict has to be overwritten.

[brownbaerchen]

Actually, in RK you don't need to compute the residual at all, so ideally this is never even called. The predict function you also don't have to overload, I think, because there is no need to fill level.f as it is not used until computed. In SDC you have this term $(Q-Q\Delta)f(u^k)$ on the right hand side, but in DIRK $Q=Q\Delta$, so to speak, and this term cancels completely. Only $Q_\Delta f(u^{k+1})$ remains, which has yet to be computed at the time of predict. So, afaik, you only need to do update_nodes for DAEs. Then, we need to figure out some inheritance magic so that you can use the existing RK methods. But that won't be too difficult. Worst case, you make a new class which inherits from the DAE RK sweeper and just references the Butcher tableaux of the existing regular RK methods.

[lisawim]

So, the original way in Runge_Kutta.py is, that RungeKutta implements the structure of the method to solve ODEs, and then each method like ForwardEuler, BackwardEuler and so on having different nodes, weights, matrix inherits from RungeKutta. Assume we have a class RungeKuttaDAE inheriting from RungeKutta, that only implements the update_nodes function. Then I see that each class defining its own nodes, weights and nodes have to be written separately for the DAE case, because each of these classes inheriting from RungeKutta.

For instance, for ForwardEuler we have that

ForwardEuler(RungeKutta)

but for DAEs we would then need

ForwardEulerDAE(RungeKuttaDAE)

but both classes would have the same attributes for the nodes, weights, and the matrix. Or did I miss something?

[brownbaerchen]

I think you could do

class ForwardEulerDAE(ForwardEuler, RungeKuttaDAE):
    pass

with RungeKuttaDAEs only implementing the update_nodes function and the nodes etc. being inherited from the regular ForwardEuler, where they are classmethods. Maybe the order of inheritance needs to be switched, I don't know. Multiple inheritance is a bit confusing, but it can be done ;)

[brownbaerchen]

On the other hand, it may be more convenient to have a single RungeKuttaDAE class that takes a regular RK method as an argument during instantiation. Than you could do

class RungeKuttaDAE(RungeKutta):
    def __init__(self, params, RKM):
        super().__init__(params)
        self.coll = RKM.get_Butcher_tableau()

    def update_nodes():
          ...

ForwardEulerDAE = RungeKuttaDAE({}, ForwardEuler)

Or something of sorts.

[lisawim]

Ah, this is what you meant by inheritance magic. :ok_hand: Thanks for discussion, I'll try that!

[lisawim]

In update_nodes you compute the values at the intermediate stages $t_{n-1} + c_i \Delta t$ by solving the implicit system $$km = f(t{n-1} + \Delta t cm, u{n-1} + \Delta t \sum^M{j=1} a{mj} k_j)$$ to get the $k_m$. So far, so good. In self.Qmat[-1, 1:-1] = weights the values of the weights $b_1,..,bM$ can be found when I see this correctly. So, for computing the value at the "final" time $t{n}$, which is the end of the time step, one has to compute the value at this point by summing up all the values at the computed intermediate stages using these weights, i.e., $$un = u{n-1} + \Delta t \sum_{j=1}^M b_j k_j$$ But in the last loop, i.e., for m = M-1 also an implicit system is solved by using the weights at the diagonal. Why? Because you don't have to solve any system there? Or did I miss something?

[brownbaerchen]

This is just out of laziness, thb. The diagonal element in the last row is always zero. You can use solve_system to get the correct solution nonetheless, even though other implementations may be more efficient.

[lisawim]

Ok, got it. :D

[lisawim]

One more question for @brownbaerchen: I already implemented BackwardEulerDAE and I'm wondering the whole time, why the solution seems to be not correct. Going several times to my code leads to no success in finding the bug and I was also doing some computations by hand to see what's going on there.. I guess I find it..

Expressing the backward Euler as Runge-Kutta method, we would get for an ODE $$k1 = f(t{n-1} + c1 \Delta t, u{n-1} + \Delta t a_{11} k_1)$$ for $c1 = 1$, $a{11}=1$ which indicates that this stage is solved implicitly (since we are looking for value $k1$ at time $t{n}$). BackwardEuler class provides

class BackwardEuler(RungeKutta):
    nodes = np.array([0.0])  # c_m
    ...

Now the update_nodes method in RungeKutta solves this stage by doing

lvl.u[m + 1][:] = prob.solve_system(
    rhs, lvl.dt * self.QI[m + 1, m + 1], lvl.u[m], lvl.time + lvl.dt * self.coll.nodes[m]
)

by having

[lvl.time + lvl.dt * self.coll.nodes[m] for m in range(M)] = [0., 0.]

Since Runge-Kutta methods are not iterative, they solve the system to the next stage $t_{n-1}+c1\Delta t$ which is $t{n-1}+\Delta t$ for backward Euler. So, my question: is it possible that BackwardEuler has to provide nodes = np.array([1.0]) instead and solves the system to a time lvl.time + lvl.dt * self.coll.nodes[m + 1] (because first node is zero anyway as stated here? For problems not using the time in their system this does not make any difference, but problems that uses the current time, it does indeed.

[brownbaerchen]

This is indeed a bug. Thanks for noticing! Just change it to 1 :)

[lisawim]

Great! No problem!