Custom solver

[KulaginVladimir]

Proposed changes

Following #673, this PR introduces:

• the possibility for the users to set the Newton solver preconditioner via F.Settings(preconditioner="...")
• the possibilty for the users to set a custom Newton solver after the model (HeatTransferProblem) initialisation.

The expected usage of a custom Newton solver is:

my_model.initialise()

class CustomSolver(fenics.NewtonSolver):
     .....

custom_solver = CustomSolver()
custom_solver.parameters["..."] = ...

my_model.h_transport_problem.newton_solver = custom_solver

my_model.run()

The same can be done for for the F.ExtrinsicTrap class:

my_model.traps = [F.ExtrinsicTrap(...), ...]
for trap in my_model.traps:
    if isinstance(trap, F.ExtrinsicTrap)
        trap.newton_solver = CustomSolver()

and for the HeatTransferProblem (but before my_model.initialise() for the steady-state case!):

my_model.T = F.HeatTransferProblem(...)
my_model.T.newton_solver = CustomSolver()

Types of changes

What types of changes does your code introduce to FESTIM?

• Bugfix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (fix or feature that would cause existing functionality to not work as expected)
• Code refactoring
• Documentation Update (if none of the other choices apply)
• New tests

Checklist

• Black formatted
• Unit tests pass locally with my changes
• I have added tests that prove my fix is effective or that my feature works
• I have added necessary documentation (if appropriate)

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 99.53%. Comparing base (4573636) to head (b9fd3eb).
Report is 190 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #734      +/-   ##
==========================================
+ Coverage   99.52%   99.53%   +0.01%     
==========================================
  Files          58       59       +1     
  Lines        2501     2577      +76     
==========================================
+ Hits         2489     2565      +76     
  Misses         12       12              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[KulaginVladimir]

One thing to mention is that the info messages, shown when log_level=20, differ with f.NewtonSolver, since f.NonlinearVariationalSolver prints additional messages. A possible workaround is to invoke f.info before sovling the problem in FESTIM.

[RemDelaporteMathurin]

@KulaginVladimir thanks for this I think it is much needed to give more flexibility.

Would you mind listing the differences - if any - from the user point of view when running the default behaviour?
Are there any differences in absolute errors, usage, etc.

[RemDelaporteMathurin]

Also tagging @Allentro as he was interested in this feature

[KulaginVladimir]

As I tested, there shouldn't appear significant differences (in errors or usage) with these changes if, for example, my_model.h_transport_problem.newton_solver is not provided manually.

During the model initialisation, a default fenics.NewtonSolver object is created and the solver parameters are set from my_model.settings: absolute_tolerance, relative_tolerance, maximum_iterations, linear_solver, and preconditioner, which were kept unchanged (the default preconditioner is set to default). Together with the introduced festim.Problem class, new changes should reproduce the solving procedure of f.NonlinearVariationalProblem and f.NonlinearVariationalSolver, used in the current version of FESTIM.

As I mentioned above, there is a change in fenics logs (i.e. when log_level=20). When NewtonSolver.solve is called, it doesn't print "Solving nonlinear variational problem.", whereas f.NonlinearVariationalSolver does.

@RemDelaporteMathurin is this what you asked? I also wonder about an example of using a custom solver. Shall we add it somewhere?

[KulaginVladimir]

So, I added fenics.begin()/fenics.end() (see decription) to increase the indentation level of info messages printed by the Newton solver.

[RemDelaporteMathurin]

I would like to keep this open while we test it thoroughly. Playing with the solver can be a bit dangerous so I want to make sure we're not merging any breaking stuff.

[RemDelaporteMathurin]

Suggested change

	options can be veiwed by print(list_krylov_solver_preconditioners()).
	options can be viewed by print(list_krylov_solver_preconditioners()).

[KulaginVladimir]

Corrected

[RemDelaporteMathurin]

Hi @KulaginVladimir my deepest apologies for taking so much time reviewing this.

I tested it and it seems to work perfectly well!

Several points:

• If users specify a custom solver for say the H transport problem, then the parameters given in Settings (absolute, relative tolerances, etc.) will be ignored. This is ok, but we should add a print statement to inform the users
• Can be done in a separate PR but it would be nice to add some documentation and maybe a short demo in the workshop?

[RemDelaporteMathurin]

Thank you so much for this awesome contribution @KulaginVladimir !!