4 p-Laplacian

Tutorial 4: p-Laplacian

In this tutorial, we will learn

Problem statement

The goal of this tutorial is to solve a nonlinear PDE in Gridap. For the sake of simplicity, we consider the $p$-Laplacian equation as the model problem. Specifically, the PDE we want to solve is: find the scalar-field $u$ such that

\[\left\lbrace \begin{aligned} -\nabla \cdot \left( |\nabla u|^{p-2} \ \nabla u \right) = f\ &\text{in}\ \Omega,\\ u = 0 \ &\text{on} \ \Gamma_0,\\ u = g \ &\text{on} \ \Gamma_g,\\ \left( |\nabla u|^{p-2}\ \nabla u \right)\cdot n = 0 \ &\text{on} \ \Gamma_{\rm N}, \end{aligned} \right.\]

with $p>2$. The computational domain $\Omega$ is the one depicted in next figure, which is the same as in the first tutorial. However, we slightly change the boundary conditions here. We impose homogeneous Dirichlet and homogeneous Neumann boundary conditions on $\Gamma_0$ and $\Gamma_{\rm N}$ respectively, and in-homogeneous Dirichlet conditions on $\Gamma_g$. The Dirichlet boundaries $\Gamma_0$ and $\Gamma_g$ are defined as the closure of the green and blue surfaces in next figure respectively, whereas the Neumann boundary is $\Gamma_{\rm N}\doteq\partial\Omega \setminus (\Gamma_0\cup\Gamma_g)$. In this example, we consider the values $p=3$, $f=1$, and $g=2$.

Numerical scheme

We discretize the problem with conforming Lagrangian FE spaces. For this formulation, the nonlinear weak form reads: find $u\in U_g$ such that $[r(u)](v) = 0$ for all $v\in V_0$. As in previous tutorials, the space $U_g$ is the set of functions in $H^1(\Omega)$ that fulfill the Dirichlet boundary conditions, whereas $V_0$ is composed by functions in $H^1(\Omega)$ that vanish at the Dirichlet boundary. The weak residual $r(u)$ evaluated at a function $u\in U_g$ is the linear form defined as

\[[r(u)](v) \doteq \int_\Omega \nabla v \cdot \left( |\nabla u|^{p-2}\ \nabla u \right) \ {\rm d}\Omega - \int_\Omega v\ f \ {\rm d}\Omega.\]

In order to solve this nonlinear weak equation, we consider a Newton-Raphson method, which is associated with a linearization of the problem in an arbitrary direction $\delta u\in V_0$, namely $[r(u+\delta u)](v)\approx [r(u)](v) + [j(u)](v,\delta u)$. In previous formula, $j(u)$ is the Jacobian evaluated at $u\in U_g$, which is the bilinear form

\[[j(u)](v,\delta u) = \int_\Omega \nabla v \cdot \left( |\nabla u|^{p-2}\ \nabla \delta u \right) \ {\rm d}\Omega + (p-2) \int_\Omega \nabla v \cdot \left( |\nabla u|^{p-4} (\nabla u \cdot \nabla \delta u) \nabla u \right) \ {\rm d}\Omega.\]

Note that the solution of this nonlinear PDE with a Newton-Raphson method, will require to discretize both the residual $r$ and the Jacobian $j$. In Gridap, this is done by following an approach similar to the one already shown in previous tutorials for discretizing the bilinear and linear forms associated with a linear FE problem. The specific details are discussed now.

Discrete model

As in previous tutorials, the first step to solve the PDE is to load a discretization of the computational domain. In this case, we load the model from the same file as in the first tutorial

using Gridap
model = DiscreteModelFromFile("../models/model.json")

As stated before, we want to impose Dirichlet boundary conditions on $\Gamma_0$ and $\Gamma_g$, but none of these boundaries is identified in the model. E.g., you can easily see by writing the model in vtk format

writevtk(model,"model")

and by opening the file "model_0" in Paraview that the boundary identified as "sides" only includes the vertices in the interior of $\Gamma_0$, but here we want to impose Dirichlet boundary conditions in the closure of $\Gamma_0$, i.e., also on the vertices on the contour of $\Gamma_0$. Fortunately, the objects on the contour of $\Gamma_0$ are identified with the tag "sides_c" (see next figure). Thus, the Dirichlet boundary $\Gamma_0$ can be build as the union of the objects identified as "sides" and "sides_c".

Gridap provides a convenient way to create new object identifiers (referred to as "tags") from existing ones. First, we need to extract from the model, the object that holds the information about the boundary identifiers (referred to as FaceLabels):

labels = FaceLabels(model)

Then, we can add new identifiers (aka "tags") to it. In the next line, we create a new tag called "diri0" as the union of the objects identified as "sides" and "sides_c", which is precisely what we need to represent the closure of the Dirichlet boundary $\Gamma_0$.

add_tag_from_tags!(labels,"diri0",["sides", "sides_c"])

We follow the same approach to build a new identifier for the closure of the Dirichlet boundary $\Gamma_g$. In this case, the boundary is expressed as the union of the objects identified with the tags "circle", "circle_c", "triangle", "triangle_c", "square", "square_c". Thus, we create a new tag for $\Gamma_g$, called "dirig" simply as follows:

add_tag_from_tags!(labels,"dirig",
  ["circle","circle_c", "triangle", "triangle_c", "square", "square_c"])

FE Space

Now, we can build the FE space by using the newly defined boundary tags.

V = FESpace(
  reffe=:Lagrangian, order=1, valuetype=Float64,
  conformity=:H1, model=model, labels=labels,
  diritags=["diri0", "dirig"])

The construction of this space is essentially the same as in the first tutorial (we build a continuous scalar-valued Lagrangian interpolation of first order). However, we also pass here the labels object (that contains the newly created boundary tags). From this FE space, we define the test and trial FE spaces

g = 1.0
V0 = TestFESpace(V)
Ug = TrialFESpace(V,[0.0,g])

Nonlinear FE problem

At this point, we are ready to build the nonlinear FE problem. To this end, we need to define the weak residual and also its corresponding Jacobian. This is done following a similar procedure to the one considered in previous tutorials to define the bilinear and linear forms associated with linear FE problems. In this case, instead of an AffineFETerm (which is for linear problems), we use a NonLinearFETerm. An instance of NonLinearFETerm is constructed by providing the integrands of the weak residual and its Jacobian (in a similar way an AffineFETerm is constructed from the integrands of the bilinear and linear forms).

On the one hand, the integrand of the weak residual is build as follows

using LinearAlgebra: norm
const p = 3
@law flux(x,∇u) = norm(∇u)^(p-2) * ∇u
f(x) = 1.0
res(u,v) = inner( ∇(v), flux(∇(u)) ) - inner(v,f)

Function res is the one representing the integrand of the weak residual $[r(u)](v)$. The first argument of function res stands for the function $u\in U_g$, where the residual is evaluated, and the second argument stands for a generic test function $v\in V_0$. Note that we have used the macro @law to construct the "constitutive law" that relates the nonlinear flux with the gradient of the solution.

On the other hand, we implement a function jac representing the integrand of the Jacobian

@law dflux(x,∇du,∇u) =
  (p-2)*norm(∇u)^(p-4)*inner(∇u,∇du)*∇u + norm(∇u)^(p-2) * ∇du
jac(u,v,du) = inner(  ∇(v) , dflux(∇(du),∇(u)) )

The first argument of function jac stands for function $u\in U_g$, where the Jacobian is evaluated. The second argument is a test function $v\in V_0$, and the third argument represents an arbitrary direction $\delta u \in V_0$. Note that we have also used the macro @law to define the linearization of the nonlinear flux.

With these functions, we build the NonLinearFETerm as follows:

trian = Triangulation(model)
quad = CellQuadrature(trian,degree=2)
t_Ω = NonLinearFETerm(res,jac,trian,quad)

We build the NonLinearFETerm by passing in the first and second arguments the functions that represent the integrands of the residual and Jacobian respectively. The other two arguments, are the triangulation and quadrature used to perform the integrals numerically. From this NonLinearFETerm object, we finally construct the nonlinear FE problem

op = NonLinearFEOperator(V,Ug,t_Ω)

Here, we have constructed an instance of NonLinearFEOperator, which is the type that represents a general nonlinear FE problem in Gridap. The constructor takes the test and trial spaces, and the FETerms objects describing the corresponding weak form (in this case only a single term).

Nonlinear solver phase

We have already built the nonlinear FE problem. Now, the remaining step is to solve it. In Gridap, nonlinear (and also linear) FE problems can be solved with instances of the type NonLinearFESolver. The type NonLinearFESolver is a concrete implementation of the abstract type FESolver particularly designed for nonlinear problems (in contrast to the concrete type LinearFESolver which is for the linear case).

We construct an instance of NonLinearFESolver as follows:

using LineSearches: BackTracking
nls = NLSolver(
  show_trace=true, method=:newton, linesearch=BackTracking())
solver = NonLinearFESolver(nls)

Note that the NLSolver function used above internally calls the nlsolve function of the NLsolve package with the provided key-word arguments. Thus, one can use any of the nonlinear methods available via the function nlsolve to solve the nonlinear FE problem. Here, we have selected a Newton-Raphson method with a back-tracking line-search from the LineSearches package.

We are finally in place to solve the nonlinear FE problem. The initial guess is a FEFunction, which we build from a vector of random (free) nodal values:

import Random
Random.seed!(1234)
x = rand(Float64,num_free_dofs(Ug))
uh0 = FEFunction(Ug,x)
uh, = solve!(uh0,solver,op)

We finish this tutorial by writing the computed solution for visualization (see next figure).

writevtk(trian,"results",cellfields=["uh"=>uh])

This page was generated using Literate.jl.