How to implement a SDC sweeper for semi-explicit DAEs?

[lisawim]

Consider a semi-explicit DAE of the form $$u'=f(u, p, t),$$ $$0 = g(u,p,t),$$ then it is reasonable that only the differential variable $u$ should be integrated, since the DAE does not include any derivative for the algebraic variables $p$. For SDC, the collocation problem to be solved for this type of DAE is then given by $$\bf{U} = f(\bf{u}_0 + \Delta t (\bf{Q}\otimes\bf{I}) \bf{U}, \bf{z}, \bf{\tau}),$$ $$0 = g(\bf{u}_0 + \Delta t (\bf{Q}\otimes\bf{I}) \bf{U}, \bf{z}, \bf{\tau}).$$ Let us call this formulation the semi-explicit collocation problem.

In pySDC, DAE problem classes uses the formulation of a general DAE of the form $$f(u', u, t)=0,$$ the method to evaluate the right-hand side of the DAE has the form eval_f(u, du, t). In this case, no distinction is made between differential variables and algebraic variables here. Using the fully-implicit sweeper fully_implicit_DAE it does not make any difference how to order the equations here. It does make sense to adopt eval_f to all other DAE problems by having different sweepers available for solving that. But using a semi-explicit sweeper to solve the semi-explicit collocation problem it indeed does make a difference how the equations are ordered by using the eval_f method. The new sweeper has to taken into account this structure.

A first idea was introducing a new attribute self.diff_nvars to all implemented semi-explicit DAEs that indicates the number of differential variables and hence the number of equations of the form $u'=f(u,p,t)$. The number of algebraic variables $p$ and thus the number of equations of the form $0=g(u,p,t)$ is then given by self.nvars - self.diff_nvars, where self.nvars is the number of unknowns of the whole system. (First version of the sweeper can be found here.) In L.u[m][:P.diff_nvars] the differential variables are stored, and L.u[m][P.diff_nvars] includes the algebraic variables for each collocation node $\tau_m$. Implementing new classes that are intended to be solved using the SemiExplicitDAE sweeper should taken the ordering of the equations taken into account, but this seems quite a bit cumbersome(?).

A small research led to DifferentialEquations.jl, where semi-explicit DAEs are also implemented in this order mentioned above, but also provides the possibility to indicate where the differential equations in the system can be found by setting

diff_nvars = [True, True, False]

Could this also be a possibility to use this for the new sweeper? However, this can result in a very cumbersome implementation for the sweeper. What would represent best practice here? Are there any other ideas to realise this? Maybe the two @pancetta and @brownbaerchen have some ideas for that issue.

[brownbaerchen]

What about a new datatype? For IMEX SDC, the datatype for f contains two meshes for the implicit and explicit part. Could you make a datatype for the solution that has two meshes, one for the differential equations and one for the algebraic ones? This could lead to more problems because many of the implemented stuff requires adding u and so on. But I guess such functions could be overloaded to act on both meshes simultaneously. In this case, you could also just evaluate f at u. So the datatype would be along the lines of

class DAE_mesh(dtype):
    def __init__(...):
        self.diff = mesh(...)
        self.alg = mesh(...)

    def __add__(self, other):
         res = DAE_mesh(...)
         res.diff = self.diff + other.diff
         res.alg = self.alg + other.alg

...

And then in the problem class

class DAE_problem(ptype):
...
    def eval_f(self, u, t):
         diff_vars = u.diff
         alg_vars = u.alg
         ...

I guess, you want level.f and level.u[:].diff to coincide, which you would have to take care of in the sweeper.

I have no idea if this is a sensible approach. I learned the other day that there are many subtle differences when doing the DAE stuff that I have no clue about. What do you think?

[pancetta]

I feel that a new datatype might be the easiest solution. We had a similar question for the particles datatype, where all of the components could have been part of a single array instead. This makes user's life miserable and is prone to errors, though. Better split it by using a specialized datatype, such as we did for particles or, as @brownbaerchen mentioned, the imex schemes.