Epsilon-embedding for differential equations - Efficient implementation

[lisawim]

The $\varepsilon$-embedding method is a tool that can be used to get insights into the behavior of numerical methods when applied to stiff ODEs and DAEs. Assume we have a stiff system of ODEs of the form $$y'=f(y, z),$$ $$\varepsilon z' = g(y,z),$$ for $0 < \varepsilon \ll 1$. The case $\varepsilon=0$ leads to the semi-explicit DAE $$y'=f(y, z),$$ $$0 = g(y,z).$$ Here, the $\varepsilon$-embedding method can propose a scheme to solve the DAE.

First, SDC is applied to the stiff ODE leading to the scheme
$$\mathbf{y}^{k+1} = \mathbf{y}_0 + \Delta t (\mathbf{Q}-\mathbf{Q}_d)\otimes\mathbf{I}_d f(\mathbf{y}^k,\mathbf{z}^k) + \Delta t \mathbf{Q}_d\otimes\mathbf{I}_d f(\mathbf{y}^{k+1},\mathbf{z}^{k+1}),$$ $$\varepsilon\mathbf{z}^{k+1} = \varepsilon\mathbf{z}_0 + \Delta t (\mathbf{Q}-\mathbf{Q}_d)\otimes\mathbf{I}_a g(\mathbf{y}^k,\mathbf{z}^k) + \Delta t \mathbf{Q}_d\otimes\mathbf{I}_a g(\mathbf{y}^{k+1},\mathbf{z}^{k+1})$$ for our well-known preconditioner $\mathbf{Q}_d$ and identity matrices $\mathbf{I}_d$ and $\mathbf{I}_a$, where $d$ is the number of differential equations of $y$ and $n_a$ is the number of differential equations of $z$ (or the number of algebraic equations in the DAE case) (GitHub does not like the $\Delta$ in the index and more than double indices..)

Setting now $\varepsilon=0$ in the SDC scheme above and making sure that the solution lie on the manifold, i.e., the numerical solution satisfies $g(\mathbf{y}^k,\mathbf{z}^k)=\mathbf{0}$ the proposed SDC scheme for the DAE case then reads
$$\mathbf{y}^{k+1} = \mathbf{y}_0 + \Delta t (\mathbf{Q}-\mathbf{Q}_d)\otimes\mathbf{I}_d f(\mathbf{y}^k,\mathbf{z}^k) + \Delta t \mathbf{Q}_d\otimes\mathbf{I}_d f(\mathbf{y}^{k+1},\mathbf{z}^{k+1}),$$ $$\mathbf{0} = g(\mathbf{y}^{k+1},\mathbf{z}^{k+1})$$ And this is very cool. Why? Because you see here, that we also have $SDC\cdot f(..)$ instead of having $f(SDC)$ in the DAE case! This means we could probably use an ODE solver like generic_implicit to treat both cases!

Phew, we've now done the theory! So, now I have a question for the implementation. I know for the first time I can make it as messy as possible on my local branch, but anyway I'll ask this question when the time will come for a corresponding PR. So, let's start..

I had a few ideas to come up with that:

1. Using only one sweeper: Since we'd like to consider both cases (stiff ODE and DAE) it might make sense to treat both in one problem class and using only the generic_implicit sweeper. Here, the DAEMesh might not make much sense (because we have an ODE and a DAE and there are only diff and alg components in the DAE case), so for the overall problem class mesh would then be used. When the SDC schemes are compared, in the DAE case rhs for the algebraic equations is zero. So, this would make generic_implicit very inefficient in terms of unnecessary computations.
2. Implementing a new sweeper: To avoid unnecessary evaluations a new sweeper could be implemented. This would probably require also the implementation of a separate DAE problem class using the DAEMesh so that the new sweeper can handle the integration only of the diff components appropriately.

Both options do not feel good. The problem classes doing in principle the same (only two different cases $\varepsilon > 0$ and $\varepsilon = 0$ here), but we would then end up with unnecessary computations. It also therefore does not make sense to implement more and more stuff and to see later, that we can save many coding lines again as we already see in the case of the MultiComponentMesh to have a core class where we can implement meshes containing components very easily, for instance.

Feel free to contribute here, I'd be very happy about that! However, there is absolutely no rush with this problem, contributions at later time are also welcome!

[brownbaerchen]

That sounds cool! It would be nice to use the existing SDC stuff to solve DAEs! I don't quite get what you need, though.

For $\epsilon=0$ you have in the first line regular SDC applied to some system of ODEs and in the second line the algebraic constraints with no apparent SDC. So my question is: Where do the algebraic constraints fit into the code?
I guess in solve_system in the problem class the algebraic constraints should be reflected somehow? Is it possible to update $y$ according to the SDC step and then invert the algebraic constraint to get $z$ inside solve_system? In this case, using existing sweepers should work "out of the box", I think. Otherwise how do you aim to apply SDC to the algebraic constraints?
As you can see, I still don't understand much about DAEs. I believe in general you cannot invert the algebraic constraints like this..?

Is your idea to make a large system of variables $x = (y, z)$, with right hand side function $F = (f, g)$ and then you solve $(1-\Delta t f)y = \mathrm{rhs}$ and $g(z) = 0$ in solve_system? Is this now invertible?

Side note: You can just make a new MultiComponentMesh with other components if it would help somehow. As long as u and f have the same mesh-derived datatype, you can just use the generic_implicit datatype, no matter the components.

[pancetta]

I may have ignored that part, but the idea would be to divide by $\epsilon$ and then just take the usual SDC sweeper, right? Then do this again for smaller $\epsilon$ and see what happens?

[brownbaerchen]

I may have ignored that part, but the idea would be to divide by ϵ and then just take the usual SDC sweeper, right? Then do this again for smaller ϵ and see what happens?

Ah and then you "numerically take the limit $\epsilon\to 0$" by setting $\epsilon=10^{-9}$ or whatever?