Lagrange multipliers etc for Dirichlet conditions

[gdmcbain]

One of the techniques for imposing Dirichlet conditions in the finite element method is Lagrange multipliers [TODO: add references]. This has been mentioned in a few previous issues [TODO: ditto].

This issue is to add a fairly simple example (plane Laplace? like ex14?) demonstrating how to implement the technique in scikit-fem.

[gdmcbain]

The discussion in #112 included the advice that

There are at least two additional ways: Nitsche's method (see interior penalty example in docs, this I'd prefer over penalty) and using mixed formulation where Dirichlet condition is natural (i.e. mixed Poisson with H(div) elements).

with which I agree; however, I thought that it might be at least a worthwhile exercise in constructing the appropriate FacetBasis to see what penalization looks like in scikit-fem. Unsurprisingly, it's easy and works well for the simple example of ex14. ee81411 Is it worth including for completeness of the discussion of boundary conditions? I'm thinking of rewriting ex14.rst to explain the various techniques for essential conditions: three so far (condense, Lagrange multiplier ab6ea3a, and now penalization).

TODO: Nitsche (which I haven't learned yet) and Raviart–Thomas (is that the intended mixed formulation? I have implemented that previously in another code but found the performance abominable, presumably because I didn't have a good preconditioner, though I presume they are out there.)

For this, I was thinking that the examples in https://kinnala.github.io/scikit-fem-docs/ are a bit flat/linear/unstructured. Would a tree make more sense than a list? Also perhaps for the three existing examples on the Stokes equations and the pending one on Navier–Stokes #145.

[kinnala]

Is it worth including for completeness of the discussion of boundary conditions? I'm thinking of rewriting ex14.rst to explain the various techniques for essential conditions: three so far (condense, Lagrange multiplier ab6ea3a, and now penalization).

Extending ex14 is a good idea. Regarding penalty method, I'd like the example to be clear on its cons:

• Not consistent (exact solution does not satisfy the penalized weak form)
• Converges optimally only for low-order basis (if I remember correctly only p=1)
• Condition number increases asymptotically faster than for standard FEM discretizations. This makes it much harder to use iterative methods for the resulting linear system as their convergence rate depends on the condition number. This can be fixed by choosing penalty parameter cleverly but fixing it will further increase the discretization error.

TODO: Nitsche (which I haven't learned yet)

I think this paper by my former colleagues is a good introduction to Nitsche's method in this context: https://math.aalto.fi/reports/a530.pdf

For this, I was thinking that the examples in https://kinnala.github.io/scikit-fem-docs/ are a bit flat/linear/unstructured. Would a tree make more sense than a list?

I agree. If you have a proposition how to categorize them I'm happy to do it.

[gdmcbain]

Actually yes I think you're right. While setting these two up in #151 was a worthwhile learning exercise for me, they'd be misleading as examples as is. (Do some of those cons apply to the Lagrange multplier method too? As I implemented it. I noted that it lost the ability to exactly represent the solution and I think I read somewhere that its condition isn't great.)

Also, besides explaining the cons, I'd like to demonstrate them. The dictum in the art of storytelling is "show don't tell"; in mathematics I think tell and show is best.

And then of course the improved alternatives (Nitsche).

I'll put this issue and the PR #151 on hold till those things are added to complete it.

Thanks very much for the review and the reference.

[gdmcbain]

I think this paper by my former colleagues is a good introduction to Nitsche's method in this context: https://math.aalto.fi/reports/a530.pdf

Thank you, §2 is exactly what I needed.

[gdmcbain]

O. K., that's ex14 redone with Nitsche's method. With that I'll reopen this.

I'll also reopen the pull request #151 so that the implementation of Nitsche's method can be looked at, but I'll leave it marked as WIP until the above-suggested discussion of relative merits is added.

[gdmcbain]

Other questions to be answered in examples and their descriptions:

• For those methods involving parameters (penalization, Nitsche's), how do the results depend on the parameter, does this depend on the size of the mesh, and how should the value be chosen.
• The shortcomings of penalization are obvious, but what about condensation versus Nitsche? Superficially:
  • They both seem to work well for this toy problem.
  • Nitsche's method requires more assembly and top-level code.
  • Condensation results in a slightly smaller stiffness matrix.

But I assume that there's more to it than that. Nitsche's method is certainly of theoretical interest in that it leads on to the concept of mortaring for multiple mismatched meshes, but although it is applicable to Dirichlet (and Robin, following Juntunen & Stenberg 2007) conditions, is it ever to be preferred over condensation?

[gdmcbain]

Another approach to essential conditions is discussed in

• Ern, A., & Guermond, J.-L. (2004). Theory and Practice of Finite Elements. Springer, §8.4.3 ‘Keeping the essential nodes’.

It's argued that condensation can be skipped with Krylov solvers provided the initial guess satisfies the essential conditions.

Might be of interest in larger applications where iterative solution is needed anyway; I haven't tried it.

[kinnala]

• For those methods involving parameters (penalization, Nitsche's), how do the results depend on the parameter, does this depend on the size of the mesh, and how should the value be chosen.

Nitsche works well with a quite large range of parameters. Penalty not so well.

But since the Nitsche parameter has an upper bound which comes from the constant in a discrete inverse inequality, you are able to calculate an optimal value for the parameter (given the element shape) from a small matrix eigenvalue problem where the size of the matrix has the same size as the number of local basis functions. I have no reference at hand and I've only seen the details explained to me once on a blackboard a few years ago.

• The shortcomings of penalization are obvious, but what about condensation versus Nitsche? Superficially:

  • They both seem to work well for this toy problem.
  • Nitsche's method requires more assembly and top-level code.
  • Condensation results in a slightly smaller stiffness matrix.

But I assume that there's more to it than that. Nitsche's method is certainly of theoretical interest in that it leads on to the concept of mortaring for multiple mismatched meshes, but although it is applicable to Dirichlet (and Robin, following Juntunen & Stenberg 2007) conditions, is it ever to be preferred over condensation?

I'd say that, in the context of setting boundary conditions, Nitsche's method will be useful when

• the boundary data is irregular
• the basis functions (Nedelec, Argyris) or the Sobolev spaces (H(curl), H^2) are so complicated that we have hard time deciding which DOFs to actually condense
• the condensation can be troublesome in practice e.g. in distributed computing (although this is probably already solved via the use of Krylov methods as you described)

Of course, I also find that it's more obviously useful in the context of mortaring.

[gdmcbain]

Thank you very much for these insights. This is gold!

[gdmcbain]

In #591 #596, an alternative to skfem.utils.condense called enforce was proposed that modifies the matrices while preserving their size. The functions have similar type-signatures. Perhaps a third function could be developed to implement the Lagrange-multiplier approach; whereas the first two shrink and preserve the size, this would augment it, with the multipliers as additional degrees of freedom. A possible name: constrain.

[kinnala]

Technically ex40 uses Lagrange multipliers to set boundary conditions. But maybe you mean enforcing the discrete DOFs to some specific value using nodewise Lagrange multiplier?

[gdmcbain]

Oh, does it? That's interesting. Even though it uses condense. I shall have to study what ibasis.get_dofs does more closely in that context. Ah, ex40.y is the Lagrange multiplier. Yes, interesting.

Another example that (presumably) internally used Lagrange multipliers for essential boundary conditions was given on 05-14; that was in the framework of minimization #719.

But yes, as per 04-01, I'm thinking about a function skfem.utils.constrain that is an alternative to condense, enforce, and penalize. Plus maybe something Nitschean after that.

[gdmcbain]

Hardly a high priority though; just for completeness. More generally I envisage that pretty much everything in a standard textbook on FEM (Ern & Guermond, Braess, the Texas series, Brenner & Scott—those being selection from my own shelf) should be doable in scikit-fem, unless there turns out to be a particular reason not to.

(Lately I've been thinking that the method of characteristics #720, for example, might not be a good fit; as Glowinski & Pan (2019) put it: “the backward method of characteristics as done in, e.g., [52, 57, 58] via the so-called Lagrange–Galerkin methodology.… Albeit conceptually simple the practical implementation of the Lagrange–Galerkin methods requires a lot of ‘savoir faire’”)

This issue could reasonably be reclassified as a discussion, if I understand the distinction correctly.

[kinnala]

Remark: C matrix in https://github.com/kinnala/scikit-fem/blob/master/docs/examples/ex07.py is essentially Nitsche's method for Dirichlet boundary conditions.