Skip to content

Unify operator application for linear and linearized problems

Currently the grid operator has the two different methods

  • jacobian_apply(z,r) (for linear problems)
  • nonlinear_jacobian_apply(x,z,r) (for linearized problems from a nonlinear equation)

for the operator application. This requires two different on-the-fly operators

  • OnTheFlyOperator
  • NonlinearOnTheFlyOperator

which call the corresponding jacobian-apply function from above on the grid operator. The same thing holds for the overlapping case. However the implementation to have two different operator application functions is tedious and error-prone to work upon.

For example matrix-free linear solver backends (like ISTLBackend_SEQ_MatrixFree_Richardson in seqistlsolverbackend.hh) should be usable in both matrix-free linear PDE-solvers and matrix-free nonlinear PDE-solvers (such as MatrixFreeNewton in matrixfreenewton.hh). Therefore the matrix-free linear solver backends have to be fed with the type of on-the-fly operator as a template argument. But these backends are usually created by the user and hence he also has to create the suitable on-the-fly operator for his problem which we definitely do not want.

There are more involved examples like matrix-free block Jacobi, block SOR or block SSOR preconditioners for DG discretizations where this separation gets really messy.

The unification will be at first tackled in EXADUNE/dune-pdelab since the interface is much more evolved there. After that the changes will be backported to PDELAB/dune-pdelab. Steffen and myself came up with the following preliminary task list:

  • Implement a single operator-apply function called jacobian_apply(x,z,r) for the grid operator variants default/, fastdg/ and onestep/ that covers both linear and linearized problems. This method takes
    • x as the linearization point
    • z as the input for the operator application
    • r as the output result
  • Deprecate the jacobian_apply(z,r) and nonlinear_jacobian_apply(x,z,r) methods in the grid operator variants default/, fastdg/ and onestep/.
  • Implement a unified on-the-fly operator called LinearizedOnTheFlyOperator and the corresponding overlapping variant OverlappingLinearizedOnTheFlyOperator which is able to handle linear and linearized problems. It checks on the grid operator whether its problem is linear or nonlinear. If it is linear, it calls jacobian_apply(z,z,r) and else it checks if the linearization point has been set prior to calling jacobian_apply(x,z,r).
  • Deprecate OnTheFlyOperator, NonlinearOnTheFlyOperatorand the corresponding overlapping variants.
  • Add the method static bool isLinear(); for the local assembler variants default/, fastdg/ and onestep/. With this the grid operator is able to check whether it implements a linear or a nonlinear problem.
  • Each local operator has to export whether it is linear or nonlinear. Therefore we add static const bool isLinear; to the LocalOperatorDefaultFlags and default it to true since this breaks much less existing code. As all the other flags it can be overwritten in the local operator.
  • For the grid operator variants default/, take the nonlinearjacobianapplyengine.hh from fastdg/. This existing engine has to be modified slightly in order to cover the linear case, too: Check if the local operator is nonlinear and if so, also load the linearization point.
  • The call switch has to be adapted to call the correct variant of the methods jacobian_apply_volume(), jacobian_apply_skeleton(), jacobian_apply_boundary() and jacobian_apply_volume_post_skeleton().

I will come up with a first WIP merge request in the next days once !5 (merged) gets accepted.