Thanks for your answer, Leo. Actually, I did just the same with the NonlinearVariationalSolver object.

In order to avoid running into new issues caused by using a NewtonSolver object instead of NonlinearVariationalSolver, e.g. because of having to solve a problem defined as NonlinearVariationalProblem and not as NonlinearProblem, I am wondering how to set the linear solver
when working with the NonlinearVariationalSolver .

Actually, I am wondering how to explicitly choose one of the solvers which can be displayed by

list_linear_solver_methods()

for the NonlinearVariationalProblem as described in the Fenics tutorial book (pp. 117–118) for linear problems.

Also, I am not sure if the method to change the solver proposed by zimtzucker is successfull as it is very similar to the command I had tried.
Is there a way to verify the type of the linear solver used or to display the NonlinearVariationalSolver’s parameters in general?

*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PCSETUP_FAILED, residual norm ||r|| = 0.000000e+00)

Apparently, a Krylov method was used although I tried to enforce the usage of a LU solver.
Actually, I am trying to avoid the error by applying a direct solution method instead of the Krylov solver.

Is there really nobody knowing how to use a direct method for nonlinear problems in Fenics? Please, I need your help.