Discontinuous Galerkin for the Nonlinear Heat Equation

Dear community,

I am looking to solve the steady heat equation with temperature-dependent conductivity using the discontinuous Galerkin method. To this end, I have used Tutorial 2.8 [1] as my starting point. The code from this tutorial runs smoothly for a simple linear test problem. My next step was to update the code to use a Newton solver (like I would for a nonlinear problem) and apply it to the same linear test problem. Using the same recipe as has previously brought success in a CG-setting, I swap the sign of the term(s) that usually go in the linear form and add them instead to the bilinear form, a. I then apply the Newton solver to solve a=0. In particular, I used the Newton solver from Tutorial 3.7 [2], which has worked well for me in the past. Unfortunately, the Newton solver diverges to infinity within roughly 10 iterations.

Any hints at how to solve this problem would be much appreciated. My code can be found on Github [3].

[1] 2.8 Discontinuous Galerkin Methods — NGS-Py 6.2.1808 documentation
[2] 3.7 Nonlinear problems — NGS-Py 6.2.2302-87-ga5a5eff3b documentation
[3] MWE_nonlinear_DG · GitHub

Best regards,

Hi Sindre,

your MWE is working with the upcoming nightly release (boundary-skeleton terms were missing in the AssembleLinearization).

However, I have to warn you. DG and non-linear assembling is not tested too much from our side, be careful. I would recommend to switch to HDG, if possible for you.


Thanks a lot, Joachim

I confirm that the fix in the nightly update solved my problem. Unfortunately, I immediately encountered another (cf. bug report [1]), so I decided to give HDG a go, although I have not used it before.

Since I am interested in nonlinear problems, I wanted to adapt the HDG example from Tutorial 2.8 to use a Newton solver (the same one I used earlier). However, when I go about this the usual way, the Newton solver complains that my “a” matrix is singular (cf. minimal not-working example [2]). It is not clear to me why this is the case, as I have not modified “a” apart from the addition of the load term. Any guidance would be much appreciated!

(Note that I am here trying to solve the same linear problem as in Tutorial 2.8, just using a Newton solver instead of the matrix inversions used in the tutorial.)

[1] Integrals over named internal faces behaves unexpectedly for DG · Issue #63 · NGSolve/ngsolve · GitHub
[2] MWE_nonlinear_HDG.py · GitHub


the fes.FreeDofs(…) requires an argument whether we want to invert on all non-constrained dofs, or only on the coupling dofs remaining after static condensation. Change this and your MWE is working:

du.data = a.mat.Inverse(fes.FreeDofs(a.condense)) * res

See also section 1.4 on static condensation.


Thanks again for your help.

For anyone looking at this in the future, in addition to updating the line mentioned by Joachim, I also had to add the following two lines immediately below that line:
du.data += a.harmonic_extension * du
du.data += a.inner_solve * res