This page was generated from unit-3.3-nonlinear/nonlinear.ipynb.

3.3 Nonlinear problems

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

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).

[1]:
# define geometry and generate mesh
from ngsolve import *
from ngsolve.webgui import *
from netgen.occ import *
shape = Rectangle(1,1).Face()
shape.edges.Min(X).name="left"
shape.edges.Max(X).name="right"
shape.edges.Min(Y).name="bottom"
shape.edges.Max(Y).name="top"
geom = OCCGeometry(shape, dim=2)
mesh = Mesh(geom.GenerateMesh(maxh=0.3))

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)\]
[2]:
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

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\):

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

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

      • Solve \(B^i \Delta u^i = -f^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:

[3]:
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:

[4]:
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.8743634219605203
Iteration   1  <A u 1 , A u 1 >_{-1}^0.5 =  0.010378681410246811
Iteration   2  <A u 2 , A u 2 >_{-1}^0.5 =  1.122416328991341e-06
Iteration   3  <A u 3 , A u 3 >_{-1}^0.5 =  1.3180540078433946e-14
[5]:
Draw(gfu,mesh,"u", deformation = True)
[5]:
BaseWebGuiScene
[6]:
Draw(gfu_it,mesh,"u", deformation = True)
[6]:
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:

[7]:
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.

We call this Newton method as well here.

[8]:
gfu.Set((x*(1-x))**4*(y*(1-y))**4) # initial guess
Newton(a,gfu,freedofs=gfu.space.FreeDofs(),maxit=100,maxerr=1e-11,inverse="umfpack",dampfactor=1,printing=True)
Newton iteration  0
err =  1.87436342196052
Newton iteration  1
err =  0.010378681410246426
Newton iteration  2
err =  1.1224163290972609e-06
Newton iteration  3
err =  1.3421225797380163e-14
[8]:
(0, 4)

A trivial problem:

As a second simple problem, let us consider a trivial scalar problem:

\[5 u^2 = 1, \qquad u \in \mathbb{R}.\]

We chose this problem and put this in the (somewhat artificially) previous PDE framework and solve with Newton’s method in a setting that you could easily follow with pen and paper:

[9]:
V = NumberSpace(mesh)
u,v = V.TnT()
a = BilinearForm(V)
a += ( 5*u*u*v - 1 * v)*dx
gfu = GridFunction(V)
gfu.vec[:] = 1
SimpleNewtonSolve(gfu,a, callback = lambda gfu : print(f"u^k = {gfu.vec[0]}, u^k**2 = {gfu.vec[0]**2}"))
print(f"\nscalar solution, {gfu.vec[0]}, exact: {sqrt(0.2)}, error: {abs(sqrt(0.2)-gfu.vec[0])}")
u^k = 1.0, u^k**2 = 1.0
Iteration   0  u^k = 0.6, u^k**2 = 0.36
<A u 0 , A u 0 >_{-1}^0.5 =  1.264911064067352
Iteration   1  u^k = 0.4666666666666667, u^k**2 = 0.2177777777777778
<A u 1 , A u 1 >_{-1}^0.5 =  0.3265986323710903
Iteration   2  u^k = 0.4476190476190476, u^k**2 = 0.20036281179138318
<A u 2 , A u 2 >_{-1}^0.5 =  0.04114755998989124
Iteration   3  u^k = 0.4472137791286727, u^k**2 = 0.20000016424254927
<A u 3 , A u 3 >_{-1}^0.5 =  0.000857426926869178
Iteration   4  u^k = 0.4472135954999957, u^k**2 = 0.20000000000003376
<A u 4 , A u 4 >_{-1}^0.5 =  3.8832745226099997e-07
Iteration   5  u^k = 0.4472135954999579, u^k**2 = 0.19999999999999998
<A u 5 , A u 5 >_{-1}^0.5 =  7.979879233426059e-14

scalar solution, 0.4472135954999579, exact: 0.4472135954999579, error: 0.0

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:

[10]:
mesh = Mesh (geom.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:

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

Now, let’s apply the Newton:

[12]:
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).")
[13]:
SolveAndVisualize()
Iteration   0  <A u 0 , A u 0 >_{-1}^0.5 =  2.7921394712043925
Iteration   1  <A u 1 , A u 1 >_{-1}^0.5 =  0.00793718884351679
Iteration   2  <A u 2 , A u 2 >_{-1}^0.5 =  6.509217393036073e-08
Iteration   3  <A u 3 , A u 3 >_{-1}^0.5 =  8.6865415467306465e-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

[14]:
nu.Set(0.01)
SolveAndVisualize(multidim=False)
Iteration   0  <A u 0 , A u 0 >_{-1}^0.5 =  1.013206169763172
Iteration   1  <A u 1 , A u 1 >_{-1}^0.5 =  0.06800962526947384
Iteration   2  <A u 2 , A u 2 >_{-1}^0.5 =  0.005083383535942226
Iteration   3  <A u 3 , A u 3 >_{-1}^0.5 =  3.713899329146616e-05
Iteration   4  <A u 4 , A u 4 >_{-1}^0.5 =  2.3705837789527196e-09
Iteration   5  <A u 5 , A u 5 >_{-1}^0.5 =  1.0974648122706111e-16
above you see the solution after the Newton solve.
[15]:
nu.Set(0.01)
SolveAndVisualize()
Iteration   0  <A u 0 , A u 0 >_{-1}^0.5 =  0.29639340338301795
Iteration   1  <A u 1 , A u 1 >_{-1}^0.5 =  0.06800962526947457
Iteration   2  <A u 2 , A u 2 >_{-1}^0.5 =  0.00508338353594237
Iteration   3  <A u 3 , A u 3 >_{-1}^0.5 =  3.7138993291481797e-05
Iteration   4  <A u 4 , A u 4 >_{-1}^0.5 =  2.370583782860462e-09
Iteration   5  <A u 5 , A u 5 >_{-1}^0.5 =  1.0544618492721814e-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).
[16]:
nu.Set(0.001)
SolveAndVisualize()
Iteration   0  <A u 0 , A u 0 >_{-1}^0.5 =  0.12675477551168057
Iteration   1  <A u 1 , A u 1 >_{-1}^0.5 =  0.027588586331091413
Iteration   2  <A u 2 , A u 2 >_{-1}^0.5 =  0.18680337968282057
Iteration   3  <A u 3 , A u 3 >_{-1}^0.5 =  0.2614577214633411
Iteration   4  <A u 4 , A u 4 >_{-1}^0.5 =  0.071300044490707
Iteration   5  <A u 5 , A u 5 >_{-1}^0.5 =  0.7737532755397406
Iteration   6  <A u 6 , A u 6 >_{-1}^0.5 =  0.6933420823135062
Iteration   7  <A u 7 , A u 7 >_{-1}^0.5 =  0.4718358798674446
Iteration   8  <A u 8 , A u 8 >_{-1}^0.5 =  0.25182949248164205
Iteration   9  <A u 9 , A u 9 >_{-1}^0.5 =  7.2528029050003795
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:

[17]:
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.1)
Draw(gfu.components[0],mesh, vectors = {"grid_size" : 25})
Newton iteration  0
err =  0.8005189933117852
Newton iteration  1
err =  0.7204486383137405
Newton iteration  2
err =  0.576613280714419
Newton iteration  3
err =  0.40447493786721095
Newton iteration  4
err =  0.2447213349570403
Newton iteration  5
err =  0.1282784565746746
Newton iteration  6
err =  0.09272826355980643
Newton iteration  7
err =  0.048350552802471004
Newton iteration  8
err =  0.01890274427188475
Newton iteration  9
err =  0.002131194384751269
Newton iteration  10
err =  2.610708792182956e-05
Newton iteration  11
err =  5.111024409936468e-09
Newton iteration  12
err =  6.087333667143062e-16
[17]:
BaseWebGuiScene

Tasks

After these explanations, here are a few suggestions for simple play-around-tasks:

  • Take the first PDE example and set up the linearization linear and bilinear forms by hand and implement a Newton solver without exploiting the convenience functions NGSolve provides to you.

  • Combine unit 3.1 and this unit and write an implicit time integration solver for the unsteady Navier-Stokes equations (start with an implicit Euler)