Redesign of solver interfaces

Registered by Anders Logg

This is a summary of a discussion with Hans Petter Langtangen, Ola Skavhaug and Johan Hake.

The design of the interface for linear solvers is not consistent with the interface for solving variational problems.

There are (at least) three different aspects of this inconsistency:

1. For linear systems, one creates a solver/algorithm: LinearSolver (or KrylovSolver), but for variational problems, one creates a problem: VariationalProblem.

2. The linear solver takes some parameter as constructor arguments (method and preconditioner), but a variational problem uses only parameters.

3. We have a shortcut function solve(A, x, b) for linear systems but not for variational problems.

Here is a suggested plan for correcting this:

1. Replace VariationalProblem with VariationalSolver.

2. Use parameters for everything, including method and preconditioner, in LinearSolver.

3. Add a solve() shortcut also for variational problems.

These fixes will take care of the inconsistencies 1-3.

In addition, we should make the following changes:

4. Be consistent with the order of arguments. In Python, we always return the output argument:

  x = solver.solve(A, b)
  u = solver.solve(a, L)

In C++, we always place the output argument first (like we already do in the assembler):

  solver.solve(x, A, b);
  solver.solve(u, a, L);

For reuse of memory, the Python interfaces should allow an optional keyword argument:

  (x = ) solver.solve(A, b, x=x)
  (u = ) solver.solve(a, L, u=u)

5. The solve() shortcut function takes parameters as optional arguments:

  x = solve(A, b, parameters={"linear_solver": "iterative", "krylov_solver": {"relative_tolerance": 1e-6}})

Blueprint information

Anders Logg
Needs approval
Marie Rognes
Series goal:
Milestone target:
milestone icon 1.0-beta
Started by
Anders Logg
Completed by
Anders Logg

Related branches



Might also want LinearSystem and VariationalProblem as placeholders for A, b and a, L, bcs. There would then be three ways to solve a problem:

  x = solve(A, b)
  x = solver.solve(A, b)
  x = LinearSystem(A, b).solve()

Hake: Why cannot LinearSolver and VariationalSolver act as this place holder?

AL: Because they are not problems, they are solvers.

Also when a solver is created it can be initiated using some default system wide parameters. These are also used for the shorter solve functions.

GNW: Some thought should also be given to re-using preconditioners, LU (symbolic) factorisations, etc. We're quite weak on this at the moment.

AL: Sounds to me like that should be a separate blueprint.

Copied from mailing list:


> Perhaps both. LinearSystem/VariationalProblem hold data, and
> LinearSolver/VariationalSolver do something.


That's what I'm thinking too now. We should update the blueprint:

We also need shortcuts such as

1. Calling LinearSolver::solve from LinearSystem::solve

    LinearSolver solver(*this);
    return solver.solve();

2. Simple solve() function that wraps everything

  x = solve(A, b, parameters)
  u = solve(a, L, parameters)

Concerning the creation of temporaries, didn't you figure out how to
get around that in creating dolfin::refine()?

AL: VariationalProblem is now a placeholder and the algorithms are implemented in solver classes.

AL: With the new VariationProblem class and the Solver classes, this is now more or less in place. Don't think it's
worth the problem of interface changes to replace solve(A, x, b) with solve(x, A, b) etc, so I consider this blueprint as implemented.


Work Items

This blueprint contains Public information 
Everyone can see this information.