Skip to content

Interface Change for Jacobian Apply Methods in Local Operators

Status quo: The jacobian apply for nonlinear problems needs the current solution x and a vector on which to apply the jacobian. Those have the same type X. See e.g. dune/pdelab/localoperator/interface.hh.

Proposal: Make it possible that they have a different type, that means x has type X and z has type Z.

Why?

If you look at the standard case of discretizing a PDE those two vectors will have the same type. Both are local vectors that belong to the ansatz-space so the status quo makes sense. For development of matrix-free preconditioners I use the whole grid operator/local operator interface to do some solving/preconditioning steps. There is a good reason to go through this interface since it avoids duplicating the whole grid iteration and assembly code.

At the moment X will either be a LocalVector (for the standard grid operator) or a ConstAliasedVectorView (for the fast-DG grid operator). For the implementation of the matrix-free solvers I have one place where I need to solve local systems of equations for the block diagonal. It makes sense to use a LocalVector as a container for the solution of this local system, as it can be passed directly into the ISTL solvers and exactly represents what we need: A vector of local degrees of freedom.

The bad thing is: I need to pass this local vector z and the current solution x to a local operator. In the fast-DG case the two types do not match anymore since z is a LocalVector and x an ConstAliasedVectorView. If we change the interface of the jacobian apply methods this will still work, since both have the accumulate and data methods that are needed in the local operator. I have tested this and it indeed solves this problem.

Are there other solutions to this problem?

  • One could try to use ConstAliasedVectorView as a type for the local solution vector z in the fast-DG case. For this one would need to implement all methods and operators needed by ISTL in the ConstAliasedVectorView. We would have two those VectorViews, one that only has local data and one that actually stores the whole global DOF vector. It might be possible to make those two work together but it sounds a bit sketchy.

  • Changing the way the LocalVector works might be another possibility (storing a pointer to some data instead of a std::vector) but that sounds even worse.

  • I could of course work around this issue by creating a LocalVector and copying data but that would defeat the whole reason of matrix-free solvers: Avoid copying data ;)

There might of course be other possibilities I didn't see.

@christi I would be really interested in opinions. The bad thing about changing the interface is changing the interface ;), but the good thing is that a ConstAliasedVectorView will always store the whole global DOF vector and a LocalVector will always hold the local degrees of freedom. If we go for another solution we will have some sort of wired data structure that may either hold only local values or all the global values depending where it is used.