4.5. 3.3 Nonlinear problems

In this unit we turn our attention to nonlinear PDE problems and the tools that NGSolve provides to simplify it.

4.5.1. A simple scalar PDE

Let us start with a simple PDE with a nonlinearity:

\[ - \Delta u + \frac13 u^3 = 10 \text{ in } \Omega \]

on the unit square \(\Omega = (0,1)^2\).

We note that this PDE can also be formulated as a nonlinear minimization problem (cf. 3.4).

# define geometry and generate mesh
from ngsolve import *
from ngsolve.webgui import *
from netgen.occ import *

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))

In NGSolve we can solve the PDE conveniently using the linearization feature of SymbolicBFI.

The BilinearForm (which is not bilinear!) needed in the weak formulation is $\( A(u,v) = \int_{\Omega} \nabla u \nabla v + 1/3 u^3 v - 10 v ~ dx \quad ( = 0 ~ \forall~v \in H^1_0) \)$

V = H1(mesh, order=3, dirichlet=[1,2,3,4])
u,v = V.TnT()
a = BilinearForm(V)
a += (grad(u) * grad(v) + 1/3*u**3*v- 10 * v)*dx

4.5.1.1. Newton’s method

In the preparation of the formulation of Newton’s method we identify the semilinear form \(A(\cdot,\cdot)\) with the corresponding operator \(A: V \to V'\) and write (assuming sufficient regularity)

• \(A(u)\) for the linear form \([A(u)](v) = A(u,v)\) and

• the derivative of \(A\) at an evaluation point \(u\) is the bilinear form \(\delta A(u)\) with $\([\delta A(u)] (w,v) = \lim_{h \to 0} \frac{ A(u+h\cdot w,v) - A(u,v)}{h}\)$

You could compute these linear and bilinear forms manually. However, NGSolve provides functions for both operations:

• The linear form \(A(u)\) for a given (vector) u is obtained from a.Apply(vec,res) (resulting in a vector res representing the linear form)

• The bilinear form (represented by the corresponding matrix) is stored in a.mat after a call of a.AssembleLinearization(vec)

Under the hood it uses its functionality to derive CoefficientFunctions symbolically.

For example, NGSolve derives the bilinear form integrand \( \frac13 u^3 v \) w.r.t. \(u\) at \(\tilde{u}\) in direction \(w\) resulting in the integrand \(\tilde{u}^2 w v\).

This allows to form the corresponding bilinear form integrals automatically for you.

To obtain a Newton algorithm we hence only need to translate the following pseudo-code formulation of Newton’s method to NGSolve. The pseudo code is:

NewtonSolve(Pseudo code):

• Given an initial guess \(u^0\)

• loop over \(i=0,..\) until convergence:

  • Compute linearization: \(A (u^i) + \delta A(u^i) \Delta u^{i} = 0\):

    • \(r^i = A (u^i)\)

    • \(B^i = \delta A(u^i)\)

    • Solve \(B^i \Delta u^i = -r^i\)

  • Update \(u^{i+1} = u^i + \Delta u^{i}\)

  • Evaluate stopping criteria

As a stopping criteria we take \(\langle A u^i,\Delta u^i \rangle = \langle A u^i, A u^i \rangle_{(B^i)^{-1}}< \varepsilon\).

Now, here comes the same thing in NGSolve syntax:

def SimpleNewtonSolve(gfu,a,tol=1e-13,maxits=10, callback=lambda gfu: None):
    res = gfu.vec.CreateVector()
    du = gfu.vec.CreateVector()
    fes = gfu.space
    callback(gfu)
    for it in range(maxits):
        print ("Iteration {:3}  ".format(it),end="")
        a.Apply(gfu.vec, res)
        a.AssembleLinearization(gfu.vec)
        du.data = a.mat.Inverse(fes.FreeDofs()) * res
        gfu.vec.data -= du
        callback(gfu)
        #stopping criteria
        stopcritval = sqrt(abs(InnerProduct(du,res)))
        print ("<A u",it,", A u",it,">_{-1}^0.5 = ", stopcritval)
        if stopcritval < tol:
            break

Let’s apply this to the previous PDE problem:

gfu = GridFunction(V)
gfu.Set((x*(1-x))**4*(y*(1-y))**4) # initial guess
gfu_it = GridFunction(gfu.space,multidim=0)
cb = lambda gfu : gfu_it.AddMultiDimComponent(gfu.vec) # store current state
SimpleNewtonSolve(gfu, a, callback = cb)

Iteration   0  <A u 0 , A u 0 >_{-1}^0.5 =  1.874663918579182
Iteration   1  <A u 1 , A u 1 >_{-1}^0.5 =  0.010379625163238133
Iteration   2  <A u 2 , A u 2 >_{-1}^0.5 =  1.1229363547662377e-06
Iteration   3  <A u 3 , A u 3 >_{-1}^0.5 =  1.3262395256549096e-14

The iterations converge to a solution in a few steps.

Draw(gfu_it,mesh,"u", deformation = True, animate = True, autoscale = True, min = 0, max = 0.73)

BaseWebGuiScene

Use the multidim-Slider to inspect the results after the iterations.

There are also some solvers shipped with NGSolve now. The ngsolve.solvers.Newton method allows you to use static condensation and is overall more refined (e.g. with damping options). For this tutorial we will mostly stay with the simple hand-crafted SimpleNewtonSolve. Here is the alternative demonstrated once, nevertheless:

from ngsolve.solvers import *
help(Newton)

Help on function Newton in module ngsolve.nonlinearsolvers:

Newton(a, u, freedofs=None, maxit=100, maxerr=1e-11, inverse='', dirichletvalues=None, dampfactor=1, printing=True, callback=None)
    Newton's method for solving non-linear problems of the form A(u)=0.
    
    Parameters
    ----------
    a : BilinearForm
      The BilinearForm of the non-linear variational problem. It does not have to be assembled.
    
    u : GridFunction
      The GridFunction where the solution is saved. The values are used as initial guess for Newton's method.
    
    freedofs : BitArray
      The FreeDofs on which the assembled matrix is inverted. If argument is 'None' then the FreeDofs of the underlying FESpace is used.
    
    maxit : int
      Number of maximal iteration for Newton. If the maximal number is reached before the maximal error Newton might no converge and a warning is displayed.
    
    maxerr : float
      The maximal error which Newton should reach before it stops. The error is computed by the square root of the inner product of the residuum and the correction.
    
    inverse : string
      A string of the sparse direct solver which should be solved for inverting the assembled Newton matrix.
    
    dampfactor : float
      Set the damping factor for Newton's method. If dampfactor is 1 then no damping is done. If value is < 1 then the damping is done by the formula 'min(1,dampfactor*numit)' for the correction, where 'numit' denotes the Newton iteration.
    
    printing : bool
      Set if Newton's method should print informations about the actual iteration like the error. 
    
    Returns
    -------
    (int, int)
      List of two integers. The first one is 0 if Newton's method did converge, -1 otherwise. The second one gives the number of Newton iterations needed.

4.5.2. Another example: Stationary Navier-Stokes:

Next, we consider incompressible Navier-Stokes equations again. This time however stationary.

Find \(\mathbf{u} \in \mathbf{V}\), \(p \in Q\), \(\lambda \in \mathbb{R}\) so that

\[\begin{align*} \int_{\Omega} \nu \nabla \mathbf{u} : \nabla \mathbf{v} + (\mathbf{u} \cdot \nabla) \mathbf{u} \cdot \mathbf{v}& - \int_{\Omega} \operatorname{div}(\mathbf{v}) p & &= \int \mathbf{f} \cdot \mathbf{v} && \forall \mathbf{v} \in \mathbf{V}, \\ - \int_{\Omega} \operatorname{div}(\mathbf{u}) q & & + \int_{\Omega} \lambda q &= 0 && \forall q \in Q, \\ & \int_{\Omega} \mu p & &= 0 && \forall \mu \in \mathbb{R}. \end{align*}\]

The domain \(\Omega\) is still \((0,1)^2\) and we prescribe homogenuous Dirichlet bnd. conditions for the velocity, except for the top boundary where we prescribe a tangential velocity. This setup is known as “driven cavity”.

Note that we use a scalar constraint to fix the pressure level that is otherwise not controlled in the presence of pure Dirichlet conditions.

We use a higher order Taylor-Hood discretization again:

mesh = Mesh (unit_square.GenerateMesh(maxh=0.05)); 
nu = Parameter(1)

V = VectorH1(mesh,order=3,dirichlet="bottom|right|top|left")
Q = H1(mesh,order=2); 
N = NumberSpace(mesh); 
X = V*Q*N
(u,p,lam), (v,q,mu) = X.TnT()
a = BilinearForm(X)
a += (nu*InnerProduct(grad(u),grad(v))+InnerProduct(grad(u)*u,v)
      -div(u)*q-div(v)*p-lam*q-mu*p)*dx

The boundary condition:

gfu = GridFunction(X)
gfu.components[0].Set(CF((4*x*(1-x),0)),
                      definedon=mesh.Boundaries("top"))

Now, let’s apply the Newton:

def SolveAndVisualize(multidim=True):
    gfu.components[0].Set(CF((4*x*(1-x),0)),
                      definedon=mesh.Boundaries("top"))
    if multidim:
        gfu_it = GridFunction(gfu.space,multidim=0)
        cb = lambda gfu : gfu_it.AddMultiDimComponent(gfu.vec) # store current state
        SimpleNewtonSolve(gfu, a, callback = cb)
    else:
        SimpleNewtonSolve(gfu, a)
    Draw(gfu.components[0],mesh, vectors = {"grid_size" : 25})
    print("above you see the solution after the Newton solve.")
    if multidim:
        Draw(gfu_it.components[0], mesh, vectors = {"grid_size" : 25})
        print("above you can inspect the results after each iteration of the Newton solve (use multidim-slider).")

SolveAndVisualize()

Iteration   0  

<A u 0 , A u 0 >_{-1}^0.5 =  2.778827565775522
Iteration   1  

<A u 1 , A u 1 >_{-1}^0.5 =  0.007917605640137049
Iteration   2  

<A u 2 , A u 2 >_{-1}^0.5 =  6.484708878259944e-08
Iteration   3  

<A u 3 , A u 3 >_{-1}^0.5 =  8.386918266419847e-16

above you see the solution after the Newton solve.

above you can inspect the results after each iteration of the Newton solve (use multidim-slider).

The problem becomes more interesting if we decrease the viscosity \(\nu\) (nu), i.e. if we increase the Reynolds number Re. Let’s consider the previous setup for decreasing values of nu

nu.Set(0.01)
SolveAndVisualize(multidim=False)

Iteration   0  

<A u 0 , A u 0 >_{-1}^0.5 =  1.01504587980676
Iteration   1  

<A u 1 , A u 1 >_{-1}^0.5 =  0.06783652426415394
Iteration   2  

<A u 2 , A u 2 >_{-1}^0.5 =  0.005059386602444661
Iteration   3  

<A u 3 , A u 3 >_{-1}^0.5 =  3.674112448114139e-05
Iteration   4  

<A u 4 , A u 4 >_{-1}^0.5 =  2.3257465459510604e-09
Iteration   5  

<A u 5 , A u 5 >_{-1}^0.5 =  9.991134246924331e-17

above you see the solution after the Newton solve.

nu.Set(0.01)
SolveAndVisualize()

Iteration   0  

<A u 0 , A u 0 >_{-1}^0.5 =  0.29526779255990737
Iteration   1  

<A u 1 , A u 1 >_{-1}^0.5 =  0.06783652426417694
Iteration   2  

<A u 2 , A u 2 >_{-1}^0.5 =  0.005059386602440515
Iteration   3  

<A u 3 , A u 3 >_{-1}^0.5 =  3.674112448106648e-05
Iteration   4  

<A u 4 , A u 4 >_{-1}^0.5 =  2.3257465526665188e-09
Iteration   5  

<A u 5 , A u 5 >_{-1}^0.5 =  9.771563700182773e-17

above you see the solution after the Newton solve.

above you can inspect the results after each iteration of the Newton solve (use multidim-slider).

nu.Set(0.001)
SolveAndVisualize()

Iteration   0  

<A u 0 , A u 0 >_{-1}^0.5 =  0.12657659604431698
Iteration   1  

<A u 1 , A u 1 >_{-1}^0.5 =  0.029727865961693037
Iteration   2  

<A u 2 , A u 2 >_{-1}^0.5 =  0.19796367204825532
Iteration   3  

<A u 3 , A u 3 >_{-1}^0.5 =  0.24952389831003563
Iteration   4  

<A u 4 , A u 4 >_{-1}^0.5 =  0.3670737701675332
Iteration   5  

<A u 5 , A u 5 >_{-1}^0.5 =  0.5937878601771654
Iteration   6  

<A u 6 , A u 6 >_{-1}^0.5 =  0.708082016226703
Iteration   7  

<A u 7 , A u 7 >_{-1}^0.5 =  1.0046393166329242
Iteration   8  

<A u 8 , A u 8 >_{-1}^0.5 =  5.882256071268309
Iteration   9  

<A u 9 , A u 9 >_{-1}^0.5 =  15.236682674134501

above you see the solution after the Newton solve.

above you can inspect the results after each iteration of the Newton solve (use multidim-slider).

Now, viscosity is so small that Newton does not converge anymore. In this case we can fix it using a small damping parameter with the ngsolve.solvers.Newton-solver:

nu.Set(0.001)
gfu.components[0].Set(CF((4*x*(1-x),0)),definedon=mesh.Boundaries("top"))
Newton(a,gfu,maxit=20,dampfactor=0.06)
Draw(gfu.components[0],mesh, vectors = {"grid_size" : 25})

Newton iteration  0

err =  0.3890921047819176
Newton iteration  1

err =  0.3657180723874925
Newton iteration  2

err =  0.3218213791090812
Newton iteration  3

err =  0.26400206635007933
Newton iteration  4

err =  0.2014072141749928
Newton iteration  5

err =  0.1408587003208587
Newton iteration  6

err =  0.08996443055452935
Newton iteration  7

err =  0.05213847470770779
Newton iteration  8

err =  0.027105267424822752
Newton iteration  9

err =  0.01246696345934633
Newton iteration  10

err =  0.004986729231632562
Newton iteration  11

err =  0.0016954532711799732
Newton iteration  12

err =  0.00047476221662914555
Newton iteration  13

err =  0.00010444911580607448
Newton iteration  14

err =  1.671189104988423e-05
Newton iteration  15

err =  1.6711985014658884e-06
Newton iteration  16

err =  6.684829241229195e-08
Newton iteration  17

err =  2.382359663796432e-15

BaseWebGuiScene

4.5.3. Nonlinear minimization problems¶

Sometimes we have to find minimizers of functionals of the form:

Find \(u \in V\) such that

\[\begin{align*} E(u)\leq E(v) \quad \forall v \in V. \end{align*}\]

As an example consider the functional:

\[\begin{align*} E(u) = \int_{\Omega} \frac12 |\nabla u|^2 + \frac {1}{12} u^4 - 10 u \, dx. \end{align*}\]

The minimum is achived when the first variation of \(E\) at \(u\) is zero, i.e. when

\[\begin{align*} \delta E(u)[v] = 0 \quad \forall v \in V. \end{align*}\]

or in other symbols

\[\begin{align*} -\Delta u + \frac13 u^3 = 10. \end{align*}\]

The problem can be set up in the following way using the Variation keyword.

V = H1(mesh, order=4, dirichlet=".*")
u = V.TrialFunction()

a = BilinearForm(V, symmetric=True)
a += Variation ((0.5*grad(u)*grad(u) + 1/12*u**4-10*u) * dx)

Now the problem becomes equivalent to the one at the beginning of the jupiter notebook.

NGSolve provides a Newton-solver for nonlinear minimization problems as well. It is called NewtonMinimization.

from ngsolve.solvers import *

gfu = GridFunction(V)


gfu.Set((x*(1-x))**4*(y*(1-y))**4)  # initial guess
NewtonMinimization(a, gfu)

Draw(gfu, mesh, "u");

Newton iteration  0
Energy:  -2.5194856851891523e-05
err =  1.8746666146522173
Newton iteration  1
Energy:  -1.7526347043107362
err =  0.010379626788788237
Newton iteration  2
Energy:  -1.7526885765106968
err =  1.1229373271047344e-06
Newton iteration  3
Energy:  -1.7526885765113243

err =  1.3250295925195858e-14