Slow convergence of Newton's method for finite elements

Question

The application is a simple non-linear advection diffusion problem (steady state) using DGFEM. My error at each iteration is given by $$ e_{n+1} = ||\mathbf{J}^{-1}(\mathbf{u}_{n})\mathbf{F}(\mathbf{u}_{n})||_{L^{2}(\Omega)} = ||\mathbf{u}_{n+1} - \mathbf{u}_{n}||_{L^{2}(\Omega)} $$

For testing purposes, I'm using an $L^{2}$ projection of the exact solution as an initial guess, i.e. solve: $$ \sum_{K\in\mathcal{M}_{h}}\int_{K}uv = \sum_{K\in\mathcal{M}_{h}}\int_{K}u_{\text{exact}}v \implies \mathbf{M}\mathbf{u}_{0} = \mathbf{f} $$

I see perfect quadratic convergence for the steady state nonlinear diffusion problem, and for most cases of nonlinear advection diffusion. For some cases it takes a VERY large number of iterations before I see quadratic convergence. It's perplexing since the initial guess is essentially the solution I want to end up at...

1) Does the advective term mess with the convergence rate of Newton's Method? Or is this just a sign that I've computed the Jacobian part of the advective terms incorrectly? I know the quadratic rate is only locally, but seeing as how the initial guess IS the solution...

2) If I start with an initial guess of $\mathbf{u}_{0}$ using the $L^{2}$ projection as defined above, should the method converge immediately or only in the case of the exact solution being exactly represented by the DG basis? Some numerical experiments point to the latter.

Here is one simple example. The model problem is $$ \nabla\cdot(\mathbf{s}(u)u) - \nabla\cdot(\kappa(u)\nabla u) = f \;\;\;\text{in}\;\Omega\subseteq\mathbb{R}^{2}\\ u = g_{D} \;\;\;\text{on}\;\partial\Omega $$

I take the exact solution $u(x,y) = x(x-1)y(y-1)$, $\mathbf{s} = [1+u,1-u]$ and $\kappa(u) = 1+u$ on the domain $\Omega = [-1,1]^{2}$. The initial guess is taken to be $\mathbf{u}_{0}$ from the $L^{2}$ projection as described above.

For basis degree p=1, here is the result:

Refinement Level 0
Iteration: 1, Error: 4.001769e+00
Iteration: 2, Error: 1.828048e+00
Iteration: 3, Error: 6.720810e+00
Iteration: 4, Error: 5.483628e+00
Iteration: 5, Error: 6.469472e+00
Iteration: 6, Error: 2.341527e+00
Iteration: 7, Error: 4.989268e+00
Iteration: 8, Error: 4.963194e+01
Iteration: 9, Error: 2.479127e+01
Iteration: 10, Error: 1.239452e+01
Iteration: 11, Error: 6.190117e+00
Iteration: 12, Error: 3.034077e+00
Iteration: 13, Error: 1.342876e+00
Iteration: 14, Error: 8.747657e-01
Iteration: 15, Error: 5.244396e-01
Iteration: 16, Error: 1.740700e-01
Iteration: 17, Error: 5.298694e-03
Iteration: 18, Error: 1.344089e-05
Iteration: 19, Error: 2.878889e-11
Iteration: 20, Error: 8.722883e-16
Iterations for convergence: 20

Refinement Level 1
Iteration: 1, Error: 1.241453e+00
Iteration: 2, Error: 4.959870e-01
Iteration: 3, Error: 1.097905e-01
Iteration: 4, Error: 6.880627e-03
Iteration: 5, Error: 5.445979e-05
Iteration: 6, Error: 3.168481e-09
Iteration: 7, Error: 1.874668e-15
Iterations for convergence: 7

Refinement Level 2
Iteration: 1, Error: 2.627395e-01
Iteration: 2, Error: 2.936349e-02
Iteration: 3, Error: 6.427106e-04
Iteration: 4, Error: 3.024645e-07
Iteration: 5, Error: 2.557842e-13
Iterations for convergence: 5

Refinement Level 3
Iteration: 1, Error: 5.830850e-02
Iteration: 2, Error: 1.113766e-03
Iteration: 3, Error: 8.567317e-07
Iteration: 4, Error: 5.108461e-13
Iterations for convergence: 4

Refinement Level 4
Iteration: 1, Error: 1.372896e-02
Iteration: 2, Error: 6.526761e-05
Iteration: 3, Error: 2.784184e-09
Iteration: 4, Error: 1.217915e-16
Iterations for convergence: 4

Refinement Level 5
Iteration: 1, Error: 3.357581e-03
Iteration: 2, Error: 4.495462e-06
Iteration: 3, Error: 1.148118e-11
Iteration: 4, Error: 1.162268e-16
Iterations for convergence: 4

Convergence in L2-norm:
1.511590 
1.809664 
2.135559 
2.080259 
2.030937 

Convergence in H1-norm:
0.748175 
1.164069 
1.155149 
1.067994 
1.029409 

Then for $p=2$, Newton's method on the coarsest mesh takes an absurd number of iterations to converge...

Refinement Level 0
Iteration: 1, Error: 9.021838e-01
Iteration: 2, Error: 6.832202e-01
Iteration: 3, Error: 2.971411e-01
...
Iteration: 683, Error: 1.499679e+00
Iteration: 684, Error: 8.019582e-01
Iteration: 685, Error: 4.119762e-01
Iteration: 686, Error: 1.051567e-01
Iteration: 687, Error: 5.846994e-03
Iteration: 688, Error: 1.780973e-05
Iteration: 689, Error: 4.898476e-11
Iteration: 690, Error: 1.755067e-15
Iterations for convergence: 690

Refinement Level 1
Iteration: 1, Error: 8.871096e-02
Iteration: 2, Error: 8.290874e-03
Iteration: 3, Error: 8.002375e-05
Iteration: 4, Error: 4.729235e-09
Iteration: 5, Error: 4.802363e-16
Iterations for convergence: 5

Refinement Level 2
Iteration: 1, Error: 1.225342e-02
Iteration: 2, Error: 9.032665e-05
Iteration: 3, Error: 1.068461e-08
Iteration: 4, Error: 1.175537e-15
Iterations for convergence: 4

Refinement Level 3
Iteration: 1, Error: 3.305251e-03
Iteration: 2, Error: 3.496646e-06
Iteration: 3, Error: 5.546949e-12
Iterations for convergence: 3

Refinement Level 4
Iteration: 1, Error: 9.875502e-04
Iteration: 2, Error: 3.438791e-07
Iteration: 3, Error: 5.191102e-14
Iterations for convergence: 3

Refinement Level 5
Iteration: 1, Error: 2.724040e-04
Iteration: 2, Error: 2.672095e-08
Iteration: 3, Error: 3.691891e-16
Iterations for convergence: 3

Convergence in L2-norm:
4.403329 
2.867100 
2.012923 
1.775206 
1.863933 

Convergence in H1-norm:
3.108488 
2.144271 
2.073729 
2.035578 
2.017033 

And for $p=4$, convergence is in a single iteration for the given initial guess. The error is just about zero as expected.

Refinement Level 0
Iteration: 1, Error: 5.726012e-12
Iterations for convergence: 1

Refinement Level 1
Iteration: 1, Error: 1.410337e-11
Iterations for convergence: 1

Refinement Level 2
Iteration: 1, Error: 2.317826e-12
Iterations for convergence: 1

Refinement Level 3
Iteration: 1, Error: 5.361553e-12
Iterations for convergence: 1

Refinement Level 4
Iteration: 1, Error: 4.668317e-13
Iterations for convergence: 1

Refinement Level 5
Iteration: 1, Error: 2.536992e-13 
Iterations for convergence: 1

Answer 1

Your observed quadratic convergence indicates that the Jacobian is likely correct. Have you looked at the solutions for your under-resolved configurations? Galerkin optimality uses the operator norm, which contains only the symmetric part, thus the solution of the discrete system could be quite different from the projection of the exact solution. This would be very likely if you use an unstable flux, or perhaps just lack a limiter, leading to solutions with numerical artifacts. What is the maximum cell Péclet number and what numerical fluxes are you using?

Answer 2

Looks like you are seeing some oscillations in the run that you encounter large number of iterations. Maybe the Jacobian itself changed too rapidly or you have a small eigenvalue. When that happens there are a few things that did help me:

1. implement a trust radius to limit your step size to a reasonably small number; ideally the trust radius will be dynamically adjusted by the program, but it is also OK to do it manually as long as convergence problem is rare.
2. shift all eigenvalues by a small amount at the beginning of the convergence. this tend to remove oscillation but you will lose quadratic convergence, so like the first one you need to scale this shift to 0 once it is on the right track.
3. add a line search step to ensure you find the minimum along the direction predicted by Newton's method. Basically, use the direction provided by Newton's method but try to get a better step size.

For me, method 3 often worked extraordinarily well, although it does take the most effort.