When should we bail from setting solver defaults?

[JHopeCollins]

We have a set of default KSP parameters and a set of default SNES parameters, defined here. These parameters get used if no parameters are passed to solve or the variational solvers.
If an "incomplete" set of solver parameters are provided (i.e. if any(k not in user_parameters for k in default_parameters.keys())) then usually we just combine the two parameter sets, giving priority to the user's parameters (I'm ignoring what happens for Real function spaces, that's a separate thing) . This is all fine and good so far.

However, if the user provides a "pc_type", then we bail from setting defaults for "ksp_type" and "pc_factor_mat_solver_type" here. The justification for this (at least for the mat solver) is that, if a user specifies pc_type: lu, then they might expect to get the PETSc default solver, i.e. PETSc (see comment on #2142).

This causes the odd situation where, if no parameters are specified, LU+MUMPS is used, and if only "pc_type": "lu" is specified, then LU+PETSc is used. This can result in difficult to trace behaviour, especially for mixed systems where MUMPS tends to perform significantly better than PETSc (see #4018). Not many people are in the habit of regularly inspecting ksp_view. Two examples where this behaviour might pop up in "user" code would be:

1. I made a new script, and didn't set any options, so I get LU+MUMPS and everything works great. Now, I start trying out different solver options, so now I specify "pc_type": "not_lu". Still fine. However, now I want to compare, and rather than using the command line, I change the parameters in code between "pc_type": "not_lu" and "pc_type": "lu". Seems like a perfectly fine thing to do, but now my LU solver is PETSc, unlike my original script without explicit parameters.

2. I am not setting any parameters, so I am getting the default MU+LUMPS. Then, I realise that my problem is symmetric, so I provide "pc_type": "cholesky". Now, I get Choleksy+PETSc, even though I didn't say anything about changing the mat solver.

I would argue that, if we think MUMPS is a good default direct solver, then we think MUMPS is a good default direct solver and we stick to that. I suspect that the silent changing of some default parameters based on other parameters probably causes more issues/confusion than users not getting LU+PETSc if they only set "pc_type": "lu". I don't think we promise to defer to PETSc default choices in the solver parameter anywhere anyway.

If we decide to stop skipping "pc_factor_mat_solver_type" and "ksp_type" if "pc_type" is provided, then this section of set_defaults could be simplified to something like the code below. I've kept the different defaults for the matfree or matnest because otherwise specifying only "mat_type": "matfree" or "mat_type": "nest" will just crash because LU needs an assembled matrix.

default_parameters = ksp_defaults | snes_defaults
if solver_parameters.get("mat_type") in {"matfree", "nest"}:
    default_parameters["ksp_type"] = "gmres"
    default_parameters["pc_type"] = "jacobi"
return default_parameters | solver_parameters

[connorjward]

What if instead we had something arguably more pythonic:

DEFAULT_SOLVER_PARAMETERS = {
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
    ....
}

def solve(..., solver_parameters=None):
    if solver_parameters is None:
        solver_parameters = DEFAULT_SOLVER_PARAMETERS
    ...

This would be a breaking change for "matfree" and "matnest" so we could have a temporary fallback and emit some deprecation warnings for that case.

This has the benefit of being totally explicit. And arguably if a user is setting mat_type then they are also very likely to be setting solver options at the same time.