Difference between LUSolver and solve

General
1

Hi all,

I have a weak for the momentum balance equation in linear elasticity and I was just wondering what is the difference between the two solving methods in fenics?

Just using solve and specifying parameters:

res = m(avg(aold, anew, alpha_m), v_) + c(pold, avg(vold, vnew, alpha_f), v_) \
       + k(pold, avg(uold, u, alpha_f), v_)
a_form = lhs(res)
L_form = rhs(res)
solve(a_form == L_form, unew, bc_u, solver_parameters={'linear_solver': 'mumps'})

versus using LUSolver: (snippet of code from Bleyer’s Code

res = m(avg(aold, anew, alpha_m), v_) + c(pold, avg(vold, vnew, alpha_f), v_) \
       + k(pold, avg(uold, u, alpha_f), v_) 
a_form = lhs(res)
L_form = rhs(res)
# Define solver for reusing factorization
K, res = assemble_system(a_form, L_form, bc_u)
solver = LUSolver(K, "mumps")
solver.parameters["symmetric"] = True
res = assemble(L_form)
[bc.apply(res) for bc in bc_u] #bc_u.apply(res)
solver.solve(K, unew.vector(), res)

Any help would be appreciated thank you!

2

A LUSolver uses a direct method: (LU factorization) to solve linear systems of the form Ax = b.
The normal backend of the LUSolvers is PETSc, which has the following options (depending on how you have installed petsc):

  • “default”
  • “umfpack”
  • “mumps”
  • “pastix”
  • “superlu”
  • “superlu_dist”
  • “petsc”
    The solver does not reassemble the matrix K when solve is called.

In your specific code, you also use assemble_system, which applies Dirichlet conditions in a symmetric fashion by applying lifting. See for instance: https://hal.archives-ouvertes.fr/hal-03741809/ for a description of the lifting procedure.

The solve method takes in ufl-forms. This means that these has to be compiled (or searched for in cache) to get the integration kernels. Boundary conditions are applied in a non-symmetric fashion (by making an identity row in the matrix, and setting the bc directly in the rhs vector).
The solve method creates a LinearVariationalSolver, that can use either direct methods (such as those above) or iterative methods with preconditioners.

1 Like
3

Hi Dokken,

Thank you so much for your reply! Just as a quick followup, when I use both methods, they give me different solutions. The solve method seems to be able to handle nonuniform and uniform meshes whereas the LUSolver method seems to handle only uniform meshes. Do you happen to know which method would give the ‘better’ solution?

Regards
Janel

4

The LUSolver should be able to handle any linear system.

You can configure the parameters in solve to give you the same solution (up to the difference in how you set boundary conditions).

Without a minimal example illustrating the two different approaches with the two different solutions, I cannot be of any further guidance.

As long as you use direct solvers (with no preconditioning) the two solvers should give the exact same result.

If you want to solve with an iterative method, you have to use solve or equivalent classes in legacy dolfin. The method to chose for these solvers depend on the PDE you are solving.

1 Like