Linear Solver Design

[tommbendall]

Linear Solver Design

The aim of our reorganisation of Gusto's linear solvers is to allow for more flexibility and rapid development of new solver strategies. This is to be achieved through:

1. using UFL to generate the linear equation set used in the solver from each PrognosticEquationSet, making it easier to add new equation sets (as they won't require a new linear solver object).
2. for the linear solver strategy to be specified through a dictionary of solver_parameters. This would facilitate non-Quasi-Newton timesteppers from using the same solver strategy as the Quasi-Newton solvers.

Part 1 was partially achieved through #649, which implemented a new general SIQNLinearSolver object.

Part 2 is tackled with #648. This implements our current linear solvers through solver parameters.

However there are going to be several pieces of technical debt and follow-up tasks, which are captured in this discussion.

---

Summary of #648

• Rather than setting up a linear solver object, solver strategies are defined through a dictionary of solver_parameters
• The default solver_parameters used in the SemiImplicitQuasiNewton stepper are the HybridisedSolverParameters that defined in the new solver_presets.py file
• The HybridisedSolverParameters covers all currently implemented PrognosticEquationSets, through an if statement
• A new AuxiliaryPC class has been implemented (why?)
• For equations which use the Schur complement to eliminate and back-substitute a thermodynamic variable, the reduced equation solved between the elimination and substitution approach needs specifying. This is done through handwritten schur_complement routines attached the the PrognosticEquationSets. This means we have essentially moved much of the hand-coded linearisations from linear_solvers.py to the schur_complement routines of the PrognosticEquationSets, so haven't yet removed them!
• The CompressibleEulerEquation case is still not covered by the infrastructure described above, so has its own special python preconditioner (I think this is because of the non-linear pressure gradient terms and the handling of the Exner pressure?)
• A routine to update_reference_profiles has been added to the PrognosticEquationSet, as this is now equation-specific

One concern that I do have is: by auto-generating the linearised equations, and specifying the solver strategy in general ways through solver parameter dictionaries, we lose the ability to read the maths by reading the code. To address this, I think we need to do the following:

1. Write some really clear documentation about the solver strategies, with explanations about how they connect to the solver parameters
2. Find a way to nicely format the output of form.__str__() so that it gives readable maths

---

Next Steps

• Make the HybridisedSolverParameters more agnostic to the PrognosticEquationSet
• Generate the Schur complement form using FML
• Find fix to issue with not updating pythonPC when reference profiles are updated
• Allow these solver parameters to be used for NonlinearVariationalSolvers (and thus in other time discretisations)
• Documentation!

[atb1995]

Very nice summary! Few comments:

1. AuxiliaryPC allows for the passing of the handwritten Schur complement form with eliminated temperature-like variable to Firedrake's HybridizationPC
2. I think the reason we cannot cover CompressibleEulerEquation with the strategy of Schur Complement + HybridizationPC is that HybridizationPC assumes a certain structure of the equations (i.e. which terms to add/replace with trace terms). The nonlinear pressure gradient term breaks the assumption so we require a bespoke hybridized preconditioner.