AuxiliaryOperatorSNES

[JHopeCollins]

If I am trying to solve F(u)=0 but for whatever reason Newton isn't working or is slow to converge, this snes type enables specifying an auxiliary nonlinear operator G(u) that can be used as a nonlinear preconditioner.

The simplest form of nonlinear preconditioning is a for a nonlinear Richardson iteration:

$$G(u^{k+1}) = G(u^{k}) - F(u^{k})$$

where the solution uhat of F(uhat)=0 is clearly a fixed point.

There's a demo for the Allen-Cahn equation where Newton fails but using a G that lags a particular term in F succeeds. AuxiliaryOpertorSNES can also be used for continuation methods like Reynolds continuation, and many other methods.

This cherrypicks the AuxiliaryOperatorSNES from #4544 because it is ready to go but the FieldsplitSNES could do with a little more stress testing/tidying up.

Pointing at release because this is strictly new functionality. The edits to firedrake/preconditioners/base.py are just updating the docstrings to be clearer about the distinction between PCBase and SNESBase.

[dham]

minor changes pablo has requested to the demo.

[pbrubeck]

Gk1 is not G(u^{k+1}) anymore. Is this the right name? How about F_aux?

[JHopeCollins]

Fair enough I guess. Gk1 - b is also not the auxiliary form though. Maybe Gk1b? It is the auxiliary form plus a forcing term (that we happen to know we will assemble Gk - Fk into later).

[pbrubeck]

It think any meaningful name should do, as long as it is not the same one.

The in-place subtraction makes it seem that we are mutating it, but that is not what is happening under the hood, so it'd be good to get rid of it.

[JHopeCollins]

I've gone with Gk1b and added a comment for what it is.

[pbrubeck]

Maybe uold, unew would be more clear here

[JHopeCollins]

uk1 and uk matches the equations describing the iteration in the class docs and the demo so they are clearer here than old/new

[pbrubeck]

Sure, with the context of the docstrings, you could tell that uk1 refers to u_{k+1}. But a subclass will not include the docstring, so the context will get lost easily in user code. If a user cargo-culted this notation, it be hard to tell whether uk1 refers to u_{k-1} or u_{k+1}.

[danshapero]

I had to puzzle over this for a moment when I was writing the demo. I agree that a different naming scheme could be beneficial. Is there a convention established elsewhere in Firedrake?

[JHopeCollins]

Is there a convention established elsewhere in Firedrake?

Not really I don't think.

In the base class the method arguments really should stay as uk1 and uk because it corresponds directly to the description of the fixed point iteration in the docstrings of the method and the class.

If uold and unew would be clearer in the demo where you don't have the docstrings, only the form method definition, then I'm fine with that.

[pbrubeck]

Is there a convention established elsewhere in Firedrake?

While single-character names are far from ideal, in Firedrake we have well established lingo. For instance, u in the solution or TrialFunction, v is the TestFunction, F is the residual Form, V is the FunctionSpace, etc. Situations where we need to distinguish between two related variables rarely arise in public API, but I can point out a few cases.

In the (non-Irksome) time-stepping demos, we define u and uold to denote the current and previous solutions.

Throughout the multigrid codebase we use uc and uf to denote the coarse and fine solutions.

[pbrubeck]

The multigrid literature often uses u^{l} and u^{l+1}. In latex or mathjax this works well because l+1 is quite explicit. If instead we opted for names like ul and ul1 in our code, as opposed to uc anf uf, we would be decreasing the readability of the code by itself.

[pbrubeck]

The issue with uk1 is that this variable name needs further explanation, due to the conflict between u^{k-1} and u^{k+1}

[pbrubeck]

Please correct me if I'm wrong, but "lagging" the negative term here is equivalent to dropping it from G.

Also having G(u^k, u^{k+1}) = F(u^k) when u^{k+1} = u^k gives you no extra source term (b).

Maybe it'd be worth adding a comment.

[JHopeCollins]

To your second point, this is the exactly the point about the stationary point in the paragraph discussing the full expression for the preconditioned iteration.

[JHopeCollins]

Please correct me if I'm wrong, but "lagging" the negative term here is equivalent to dropping it from G.

Once you expand the full Gk1 = Gk - Fk yes that term cancels in Gk1 and Gk. But I think it is better to include it in the definition of G here because it keeps the equivalence of the fixed point iteration with G(u^{k+1}; u^{k}) = 0.

[JHopeCollins]

Also having G(u^k, u^{k+1}) = F(u^k) when u^{k+1} = u^k gives you no extra source term (b).

Yes for this specific choice of G this is true. But we cannot rely on that in general because it's perfectly acceptable for it not to be true.

[pbrubeck]

The source term b is zero for this choice of F. But one could also have obtained G by dropping the negative term in F, in which case b is nonzero, but yet you get an identical algorithm.

I think we should comment on the alternative choice of dropping the negative term.

[pbrubeck]

Also this comment would be helpful to let users know that the dependence of G on the old state is completely optional

[danshapero]

I think explaining at this level of detail will only confuse someone who's new to NPC. For that matter, I'd rather change some of the text that Josh already added to the demo back to considering only the particular case where G(u_k; u_k) = F(u_k). In that case the right-hand side of the fixed-point iteration cancels and you end up solving just G(u; u_k) = 0. This is often sufficient to obtain a good NPC and for an introduction I think that's more important to communicate than anything else. You can show people later that it isn't a necessary condition. The full details are already in the docstring for AuxOpSNES.

[pbrubeck]

Sure, I'm not too fussy about this. You might have a point, some people might get confused if we casually mention that dropping the term that depends on the old state from both sides of the equation gives an equivalent nonlinear system

[pbrubeck]

Now G(u) means a different thing. A better name would be E(u)

[danshapero]

I've changed this as you suggest. It would be a bigger change but I could instead define E(u) at the very beginning and work almost entirely from this, defining F = firedrake.derivative(E, u).

[pbrubeck]

If you do that, then the demo does not run in complex mode

[JHopeCollins]

Is Allen-Cahn ill posed in complex arithmetic then? Surely if the derivation is not well posed then it doesn't make sense to solve the equation in complex mode. We can add the skipcomplex marker to the demo in the test suite.

[danshapero]

Lev Landau has entered the chat. Is it a small change to the algebra in the form to make this work in complex mode? If that's more involved then I say we punt on complex mode for later but we do want to get to that eventually.

[pbrubeck]

It's not that the PDE is not well posed in complex arithmetic, when we symbolically differentiate inner(u,u)*dx twice to get the Jacobian we get inner(test, trial)*dx + inner(trial, test)*dx, and the form assembler requires only the test function to be conjugated.

[pbrubeck]

We could try anderson as well, it reduces iterations down to 11 for this example

Suggested change

	"snes_type": "nrichardson",
	"snes_type": "anderson",
	"snes_anderson_m": 2,

[danshapero]

I considered this but wanted to keep it simple at first. I could add Anderson or NGMRES at the end. But I was planning another demo showing more advanced techniques using the Navier-Stokes equations in a later PR.

[JHopeCollins]

How about just a sentence along the lines of "we could improve the convergence by swapping to an accelerated SNES type such as Anderson or ngmres, which you may want to experiment with in practice/in your own application"

[danshapero]

Done and I added a bit more to mention that NGMRES has more choices to make.

[pbrubeck]

This could go in the manual. This class is not going to be instantiated ever by a user.

[JHopeCollins]

But it is still a user facing class. This will get rendered in the API docs

[pbrubeck]

Regardless of what we decide to leave in the class docstring, we should be adding an entry for AuxiliaryOperatorSNES in this manual section.

I think this level of detail could live in that other section, and we could keep a brief summary in the docstring, with a cross-reference.

We documented things in a similar way when we introduced the quad_scheme kwarg for FunctionSpace in https://github.com/firedrakeproject/firedrake/pull/4607/changes

[pbrubeck]

Suggested change

	appctx=ctx.appctx, options_prefix=prefix)
	appctx=ctx.appctx,
	options_prefix=prefix,
	pre_apply_bcs=ctx.pre_apply_bcs,
	)

[pbrubeck]

In the implementation, G depends on both u^{k+1} and u^{k}. This dependence should be made evident in the docstring.

[pbrubeck]

There is one note in form that explains this, but I think the main mathematical description in this docstring should at least suggest the possibility of G depending on the old state.

[pbrubeck]

The demo refers to G as G_{k}(u) = G(u; u_{k}), should we adopt that convention here? It'd be good to unify notation.

[pbrubeck]

Now we don't test for complex by default. This will break in complex mode

Suggested change

	F = (eps * inner(grad(u), grad(v)) + (u**3 - u) * v) * dx
	F = (eps * inner(grad(u), grad(v)) + inner(u**3 - u, v)) * dx

[pbrubeck]

Same complex issue

Suggested change

	return (eps * inner(grad(u), grad(v)) + (u**3 - u_k) * v) * dx, bcs
	return (eps * inner(grad(u), grad(v)) + inner(u**3 - u_k, v)) * dx, bcs

[pbrubeck]

Just a thought. What is so special about G depending on the current state? It might as well depend on any other coefficient in the original problem. Should we instead of having the form(snes, uk, uk1, v) signature drop uk and promote grabbing the current state as uk = appctx["state"]? This aligns with the mechanism we have to access the full state for AssembledPC within a fieldsplit.