linear solver choice

Question

I've got the following piece of code

problem = LinearVariationalProblem(aLin, L, u, bc)
solver = LinearVariationalSolver(problem)
solver.parameters["linear_solver"] ="lu"
solver.solve()

and no matter what I select for the linear solver, I see in the command line:

Solving linear variational problem.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Solving linear system of size 289 x 289 (uBLAS Krylov solver).
*** Warning: Conjugate-gradient method not yet programmed for uBLASKrylovSolver. Using GMRES.
Krylov solver converged in 4 iterations

For the problem that I am interested in (studying a particular algorithm), I would like to eliminate the numerical linear algebra error, either by using a direct solver, or by performing more GMRES iterations.

Comments

you can also get FEniCS to tell you what it has installed, so that you can figure out in advance what you have to work with:

list_krylov_solver_methods()
list_krylov_solver_preconditioners()
list_linear_algebra_backends()
list_linear_solver_methods()
list_lu_solver_methods()

So for a Krylov solver, you can then set the iterations or relative tolerance like this (or, even monitor the convergence on the screen):

solver = KrylovSolver("gmres", "ml_amg")
solver.parameters["relative_tolerance"] = 5e-6
solver.parameters["maximum_iterations"] = 1000
solver.parameters["monitor_convergence"] = True

---

How do I use that in connection with the linear variational problem? I've been using the command:

problem = LinearVariationalProblem(aLin, L, u, bc)
solver = LinearVariationalSolver(problem)

but if I try to use:

solver.parameters["relative_tolerance"] = 5e-6
solver.parameters["monitor_convergence"] = True

It gives me an error that no such parameter exists.

---

In my example - I was specifically using a Krylov solver, but for your example your parameter structure has an additional layer in it. You would need to change them like this:

solver.parameters["krylov_solver"]["relative_tolerance"] = 5e-6
solver.parameters["krylov_solver"]["monitor_convergence"] = True

The parameter structure is designed to look like a dictionary in Python, if you want to examine what is in there, you can try something like this from the command line:

>>> for item in solver.parameters.items(): print item
... 
('linear_solver', 'default')
('preconditioner', 'default')
('print_matrix', False)
('print_rhs', False)
('symmetric', False)
('krylov_solver', <dolfin.cpp.common.Parameters; proxy of <Swig Object of type 'std::shared_ptr< dolfin::Parameters > *' at 0x10bef0fc0> >)
('lu_solver', <dolfin.cpp.common.Parameters; proxy of <Swig Object of type 'std::shared_ptr< dolfin::Parameters > *' at 0x10bef0f30> >)

Answer 1

You need to configure DOLFIN with a direct solver. The easiest way to do this is to configure DOLFIN with PETSc, and enable a direct solver in PETSc, e.g. using the PETSc configure options

--download-suitesparse --download-superlu --download-superlu_dist

Comments

My dolfin installation comes from macports, and it is not configured with petsc. However, it is configured with suitesparse(umfpack), and I get

>>> list_linear_solver_methods()
Solver method  |  Description                                               
----------------------------------------------------------------------------
default        |  default linear solver                                     
umfpack        |  UMFPACK (Unsymmetric MultiFrontal sparse LU factorization)
cg             |  Conjugate gradient method                                 
gmres          |  Generalized minimal residual method                       
bicgstab       |  Biconjugate gradient stabilized method                    

But if I try setting the solver to be umfpack, I still see:

Solving linear variational problem.
Solving linear system of size 81 x 81 (uBLAS Krylov solver).