A Multigrid Optimization Algorithm for the Numerical Solution of Quasilinear Variational Inequalities Involving the $p$-Laplacian

TL;DR

This paper introduces a multigrid optimization algorithm for efficiently solving quasilinear variational inequalities involving the p-Laplacian, with applications to viscoplastic fluid flow simulations.

Contribution

It develops a novel MG/OPT algorithm with a Huber regularization and finite element discretization, proving its global convergence for nonsmooth energy functionals.

Findings

The algorithm is globally convergent for the class of problems considered.

Numerical experiments demonstrate high efficiency in simulating viscoplastic fluid flows.

The approach effectively handles nonsmooth energy functionals involving the p-Laplacian.

Abstract

In this paper we propose a multigrid optimization algorithm (MG/OPT) for the numerical solution of a class of quasilinear variational inequalities of the second kind. This approach is enabled by the fact that the solution of the variational inequality is given by the minimizer of a nonsmooth energy functional, involving the $p$-Laplace operator. We propose a Huber regularization of the functional and a finite element discretization for the problem. Further, we analyze the regularity of the discretized energy functional, and we are able to prove that its Jacobian is slantly differentiable. This regularity property is useful to analyze the convergence of the MG/OPT algorithm. In fact, we demostrate that the algorithm is globally convergent by using a mean value theorem for semismooth functions. Finally, we apply the MG/OPT algorithm to the numerical simulation of the viscoplastic flow of Bingham, Casson and Herschel-Bulkley fluids in a pipe. Several experiments are carried out to show the efficiency of the proposed algorithm when solving this kind of fluid mechanics problems.

Tables

Table 1. Table 1: Results of the resolution of problem ( 38 ) with p = 1.75 𝑝 1.75 p=1.75 , g = 0.2 𝑔 0.2 g=0.2 , γ = 10 3 𝛾 superscript 10 3 \gamma=10^{3} and f = 1 𝑓 1 f=1 .

MG/OPT ( ${\nu }_1={\nu }_2=2$)					Preconditioned descent algorithm
g	It.	$Er{r}_s$	Plug Flow vel.	$Er{r}_{pf}$	It.	$Er{r}_s$	Plug Flow vel.	$Er{r}_{pf}$
	1	9.69e-09	0.169	1.42e-04	5	3.36e-09	0.1700	2.61e-05
0	2	2.44e-14	0.1700	3.57e-05	9	4.97e-14	0.1700	3.57e-05
	3	1.08e-18	0.1700	3.59e-05	12	1.0e-18	0.1700	3.59e-05
	1	4.29e-05	5.91e-02	7.5e-03	6	3.09e-05	5.94e-02	7.7e-03
	3	3.94e-06	5.36e-02	1.9e-03	31	3.99e-06	5.35e-02	1.9e-03
0.2	5	2.82e-06	5.33e-02	1.7e-03	40	2.78e-06	5.33e-02	1.6e-03
	7	1.72e-06	5.30e-02	1.3e-03	51	1.78e-06	5.29e-02	1.3e-03
	9	1.00e-06	5.27e-02	1.1e-03	61	1.06e-06	5.27e-02	1.0e-03
	1	1.6e-04	1.76e-02	1.3e-02	5	1.13e-04	1.5e-02	1.1e-02
	3	1.08e-05	6.8e-03	2.0e-03	16	1.02e-05	6.7e-03	2.7e-03
0.4	5	3.80e-06	5.7e-03	1.7e-03	35	3.73e-06	5.7e-03	1.7e-03
	7	1.86e-06	5.3e-03	1.3e-03	49	1.81e-06	5.0e-03	1.2e-03
	9	1.00e-06	5.0e-03	1.0e-03	58	1.09e-06	5.0e-03	1.0e-03

Table 2. Table 2: Results of the resolution of problem ( 38 ) with p = 2 𝑝 2 p=2 , γ = 10 3 𝛾 superscript 10 3 \gamma=10^{3} and f = 1 𝑓 1 f=1 .

MG/OPT ( ${\nu }_1={\nu }_2=2$)					Preconditioned descent algorithm
g	It.	$Er{r}_s$	Plug Flow vel.	$Er{r}_{pf}$	It.	$Er{r}_s$	Plug Flow vel.	$Er{r}_{pf}$
	1	1.0e-03	4.6e-02	3.63e-02	4	1.1e-03	5.36e-02	4.36e-02
	3	7.38e-05	1.92e-02	9.2e-03	12	6.74e-05	2.01e-02	1.01e-02
0.4	5	4.47e-05	1.70e-02	7.0e-03	15	4.08e-05	1.74e-02	7.4e-03
	7	2.96e-05	1.61e-02	6.1e-03	25	2.93e-05	1.62e-02	6.2e-03
	9	2.59e-05	1.58e-02	5.8e-03	33	2.51e-05	1.57e-02	5.7e-03

Table 3. Table 4: Results of the resolution of problem ( 38 ) with p = 5 𝑝 5 p=5 , γ = 10 3 𝛾 superscript 10 3 \gamma=10^{3} and f = 1 𝑓 1 f=1 .

MG/OPT ( ${\nu }_1={\nu }_2=2$)					Preconditioned descent algorithm
g	It.	$Er{r}_s$	Plug Flow vel.	$Er{r}_{pf}$	It.	$Er{r}_s$	Plug Flow vel.	$Er{r}_{pf}$
	1	6.67e-05	4.92e-01	1.63e-02	4	6.91e-05	4.91e-01	1.69e-02
0.1	2	6.38e-07	5.03e-01	5.6e-03	9	2.78e-07	5.01e-01	7.1e-03
	3	9.31e-09	5.02e-01	6.7e-03	51	8.98e-09	5.01e-01	7.6e-03

Table 4. Table 5: Results of the resolution of problem ( 38 ) with γ = 10 3 𝛾 superscript 10 3 \gamma=10^{3} and f = 1 𝑓 1 f=1 .

No.Grids	No. Nodes	It	Plug flow	$Er{r}_{pf}$	Rel CPU Time
1	8321	10	1.55e-02	5.03e-04	1
2	2113	2	1.51e-02	4.00e-05	0.93
3	545	3	1.51e-02	2.95e-05	0.75
4	145	4	1.51e-02	4.84e-05	0.70
5	41	5	1.50e-02	5.99e-06	0.69

Equations

$${\int }_{\Omega }|\nabla u{|}^{p-2}(\nabla u,\nabla (v-u))dx+g{\int }_{\Omega }|\nabla v|dx-g{\int }_{\Omega }|\nabla u|dx\ge {\int }_{\Omega }f(v-u)dx,\forall v\in W_0^{1,p}(\Omega ),$$

$$\underset{u\in W_0^{1,p}(\Omega )}{\min }J(u):=\frac{1}{p}{\int }_{\Omega }|\nabla u{|}^{p}dx+g{\int }_{\Omega }|\nabla u|dx-{\int }_{\Omega }fudx.$$

$$J^{\prime }(x)(h):=\underset{t\to 0^{+}}{\lim }\frac{J(x+th)-J(x)}{t}$$

$${\partial }_{B}J(x)=\left\{M\in {\mathbb{R}}^{m\times n}:\exists (x_k)\subset D_J:x_k\to x,J^{\prime }(x_k)\to M\right\}$$

$$\Vert J(x+h)-J(x)-J^{\prime }(x)(h)\Vert =o(\Vert h\Vert )\text{as}h\to 0.$$

$$J^{\prime }(x)(h)=Mh.$$

$$\underset{M\in \partial J(x+h)}{\sup }\Vert Mh-J^{\prime }(x)(h)\Vert =o(\Vert h\Vert )\text{as}h\to 0.$$

$$\underset{M\in \partial J(x+h)}{\sup }\Vert J(x+h)-J(x)-Mh\Vert =o(\Vert h\Vert )\text{as}h\to 0.$$

$$\Vert J^{\prime }(x+h)(h)-J^{\prime }(x)(h)\Vert =o(\Vert h\Vert )\text{ as }h\to 0.$$

$$\underset{h\to 0}{\lim }\frac{\Vert J(x+h)-J(x)-J^{\circ }(x+h)h\Vert }{\Vert h\Vert }=0.$$

$$\partial J(x)={\partial }_{B}J(x)=\{J^{\prime }(x)\}$$

$$\vec{z}\to \max (g,\gamma |\vec{z}|)$$

$$\chi_{A_{\gamma}}(\vec{z})=\begin{cases}\begin{array}[]{ll}1,&\text{ if }\vec{z}\in A_{\gamma},\\ 0,&\text{ if }\vec{z}\in X\setminus A_{\gamma},\end{array}\end{cases}$$

$$J(x_1)=J(x_0)+{\int }_0^1{J}^{\prime }(x_0+t(x_1-x_0))(x_1-x_0)dt$$

$${\psi }_{\gamma }:z\to {\psi }_{\gamma }(z)=\{\begin{array}{ll}g|z|-\frac{g^2}{2\gamma }& if|z|>\frac{g}{\gamma }\\ \frac{\gamma }{2}|z{|}^2& if|z|\le \frac{g}{\gamma }.\end{array}$$

$$\underset{u\in W_0^{1,p}(\Omega )}{\min }{J}_{\gamma }(u):=\frac{1}{p}{\int }_{\Omega }|\nabla u{|}^{p}dx+{\int }_{\Omega }{\psi }_{\gamma }(\nabla u)dx-{\int }_{\Omega }fudx.$$

$$V_h:=\{v_h\in C({\bar{\Omega }}_h):v_h{|}_{T_i}\in {\mathbb{P}}_1,\forall T_i\in {\Omega }_h\},$$

$$V_h^0=W_0^{1,p}(\Omega )\cap V_h$$

$$\underset{u_h\in V_h^0}{\min }{J}_{\gamma ,h}(u_h):=\frac{1}{p}{\int }_{{\Omega }_h}|\nabla u_h{|}^{p}dx+{\int }_{{\Omega }_h}{\psi }_{\gamma }(\nabla u_h)dx-{\int }_{{\Omega }_h}f{u}_{h}dx.$$

$$J_{\gamma ,h}(u_h):={\mathcal{F}}_h(u_h)+{\mathcal{G}}_{\gamma ,h}(\nabla u_h),$$

$${\mathcal{F}}_h(u_h):=\frac{1}{p}{\int }_{{\Omega }_h}|\nabla u_h{|}^{p}dx-{\int }_{{\Omega }_h}f{u}_{h}dx\text{mbox}and{\mathcal{G}}_{\gamma ,h}(\nabla u_h):={\int }_{{\Omega }_h}{\psi }_{\gamma }(\nabla u_h)dx.$$

$${\mathcal{F}}_h^{\prime }(u_h)v_h={\int }_{{\Omega }_h}|\nabla u_h{|}^{p-2}\nabla u_h\cdot \nabla v_{h}dx-{\int }_{{\Omega }_h}f{v}_{h}dx,\forall v_h\in V_h^0$$

$$\begin{array}[]{rll}\mathcal{F}^{\prime\prime}_{h}(u_{h})(v_{h},w_{h})=&\displaystyle\int_{\Omega_{h}}|\nabla u_{h}|^{p-2}\nabla v_{h}\cdot\nabla w_{h}\,dx\vspace{0.2cm}\\ &+(p-2)\int_{\Omega_{h}}|\nabla u_{2}|^{p-4}(\nabla u_{h}\cdot\nabla v_{h})(\nabla u_{h}\cdot\nabla w_{h}),\,\,\,\forall v_{h},w_{h}\in V^{0}_{h}.\end{array}$$

$$J_{\gamma ,h}^{\prime }(u_h)v_h:={\int }_{{\Omega }_h}|\nabla u_h{|}^{p-2}\nabla u_h\cdot \nabla v_{h}dx+g{\int }_{{\Omega }_h}\frac{\gamma (\nabla u_h\cdot \nabla v_h)}{\max (g,\gamma |\nabla u_h|)}dx-{\int }_{{\Omega }_h}f{v}_{h}dx,\forall v_h\in V_h^0.$$

$$\begin{array}[]{lll}(J^{\prime}_{\gamma,h})^{\circ}(u_{h})(v_{h},w_{h}):=\displaystyle\int_{\Omega_{h}}|\nabla u_{h}|^{p-2}\nabla v_{h}\cdot\nabla w_{h}\,dx\vspace{0.2cm}\\ \hskip 56.9055pt+(p-2)\displaystyle\int_{\Omega_{h}}|\nabla u_{h}|^{p-4}(\nabla u_{h}\cdot\nabla v_{h})(\nabla u_{h}\cdot\nabla w_{h})\,dx\vspace{0.2cm}\\ \hskip 85.35826pt+\displaystyle\int_{A_{\gamma}}g\frac{(\nabla v_{h}\cdot\nabla w_{h})}{|\nabla u_{h}|}\,dx-\displaystyle\int_{A_{\gamma}}g\frac{(\nabla u_{h}\cdot\nabla v_{h})(\nabla u_{h}\cdot\nabla w_{h})}{|\nabla u_{h}|^{3}}\,dx\vspace{0.2cm}\\ \hskip 156.49014pt+\displaystyle\int_{\Omega_{k-1}\setminus A_{\gamma}}\gamma(\nabla v_{h}\cdot\nabla w_{h})\,dx,\,\,\,\forall v_{h},w_{h}\in V^{0}_{h}.\end{array}$$

$${\mathcal{G}}_{\gamma ,h}^{\prime }(\nabla u_h)v_h=g{\int }_{A_{\gamma ,h}}\frac{\nabla u_h\cdot \nabla v_h}{|\nabla u_h|}dx+g{\int }_{{\Omega }_h\setminus A_{\gamma ,h}}\gamma (\nabla u_h\cdot \nabla v_h)dx,\forall v_h\in V_h^0,$$

$$A_{\gamma ,h}=\{x\in {\Omega }_h:\gamma |\nabla u_h(x)|\ge g\}.$$

$${\mathcal{G}}_{\gamma ,h}^{\prime }(\nabla u_h)v_h:=g{\int }_{{\Omega }_h}\frac{\gamma (\nabla u_h\cdot \nabla v_h)}{\max (g,\gamma |\nabla u_h|)}dx,\forall v_h\in V_h^0.$$

$$J_{\gamma ,h}^{\prime }(u_h)v_h={\mathcal{F}}_h^{\prime }(u_h)v_h+{\mathcal{G}}_{\gamma ,h}^{\prime }(\nabla u_h)v_h,$$

$$J_{\gamma ,h}^{\prime }(u_h)v_h:={\int }_{{\Omega }_h}|\nabla u_h{|}^{p-2}\nabla u_h\cdot \nabla v_{h}dx+g{\int }_{{\Omega }_h}\frac{\gamma (\nabla u_h\cdot \nabla v_h)}{\max (g,\gamma |\nabla u_h|)}dx-{\int }_{{\Omega }_h}f{v}_{h}dx\forall v_h\in V_h^0.$$

Full text

**Keywords: ** Multigrid methods. Variational inequalities of the second kind. $p$-Laplacian. Preconditioned descent algorithms. Viscoplastic Fluids.

**AMS Subject Classification: ** 65N55, 65K15, 65N30, 76A05.

1 Introduction

In this paper, we focus on the development of a multigrid algorithm for the fast finite element solution of a class of quasilinear variational inequalities of the second kind. The main idea is the application of an efficient multigrid approach to this kind of problems which, typically, leads us to the solution of large systems.

Let $\Omega$ be an open and bounded set in $\mathbb{R}^{n}$ with Lipschitz boundary $\partial\Omega$. We are concerned with the numerical solution of the following class of quasilinear variational inequalities of the second kind: find $u\in W_{0}^{1,p}(\Omega)$ such that

[TABLE]

where $1<p<\infty$, $g>0$ and $f\in L^{q}(\Omega).$ Here, $q=\frac{p}{p-1}$ stands for the conjugate exponent of $p$. It is known that these variational inequalities correspond to a first order necessary optimality condition for the following class of nonsmooth optimization problems.

[TABLE]

Consequently, we focus on the fast solution of this optimization problem. The existence and uniqueness of solutions for this problem has been analyzed and verified in previous contributions, such as [12].

The variational inequalities under study provide a versatile tool in the study of a class of free boundary problems which arise in the modelling of complex fluids and materials. In fact, diverse problems including the flow of viscoplastic materials, the flow of electro- and magneto-rheological fluids and phenomena in glaciology have been successfully simulated by this kind of models ([12, 28, 17]).

Several approaches have been proposed for the numerical solution of problems like (1). In [16] an Augmented Lagrangian method is implemented for the numerical simulation of the flow of viscoplastic materials. In [12], a preconditioned descent algorithm is proposed and analyzed both in finite and infinite dimension settings. Regarding the use of multigrid algorithms, in [18, 19] the author proposes algorithms for variational inequalities of the first and the second kind, based on extended relaxation methods. In [20], the author proposes a multigrid algorithm for variational inequalities of the second kind using a combination of convex minimization with constrained Newton linearization. However, these contributions focus on variational inequalities involving linear elliptic operators such as the Laplacian.

The multigrid approach is a very appealing way to develop fast solution algorithms for the numerical approximation of (1). In fact, the numerical solution of this kind of problems usually involves the resolution of large linear and nonlinear systems. Since these systems are computationally expensive to solve, the multigrid algorithms provide an efficient way to handle the large systems generated when discretizing the problem. Furthermore, it is natural to look for an algorithm which, in the context of the multigrid approximation, focus on the direct optimization of the energy functional.

The multigrid optimization method (MG/OPT) corresponds to a nonlinear programming adaptation of the full approximation storage (FAS) scheme. This approach is proposed, for instance in [23, 25, 31], as an effective tool for large scale optimization problems. This algorithm works with different discretization levels of the optimization problem and takes advantage of the coarse problems to generate search directions for the finer problems. Similar approaches have been used for problems involving quasilinear operators, such as the $p$-Laplace operator (see [2] and the references therein), but to the best of our knowledge, there are no contributions proposing a MG/OPT algorithm for variational inequalities of the second kind involving this kind of operators.

In this paper, we propose and analyze an MG/OPT algorithm to compute the finite element solution of a Huber regularized version of (1). Considering the structure of the optimization problem, specifically the low regularity of the functional, we use a class of preconditioned descent algorithms proposed in [12] as underlying optimization methods or smoothers. Further, the low regularity of the functional prevents us from doing a classical analysis of convergence. Therefore, we perform the convergence analysis of the MG/OPT algorithm by using a mean value theorem for Bouligand differentiable functions, which is also applicable for semismooth functions. Finally, we present a comprehensive numerical experimentation focused on the numerical simulation of viscoplastic materials. Specifically, we focus on the flow of these materials through the cross-section of a pipe.

Let us mention that, although the method developed in this article is concerned with variational inequalities of the second kind involving the nonsmooth term $\int_{\Omega}|\nabla u|\,dx$, the results can be extended to other variational inequalities of the second kind.

The paper is organized as follows: In section 2, we present several results on generalized differentiability, which will be used to analyze the convergence of the multigrid method. Since the problem is nonsmooth, in section 3 we propose a local regularization for the objective functional. Further, we present the finite element discretization of the problem. The MG/OPT method is presented in section 4, whereas the convergence of the algorithm is discussed in section 5. In section 6, a brief explanation of the underlying optimization and line search algorithms is presented. In section 7, we analyze the behaviour of the proposed methodology when applied to the numerical simulation of viscoplastic flow. We perform several experiments in order to show the main features of the algorithm. Finally, in section 8, we outline conclusions on this work and discuss future contributions.

2 Preliminaries on Generalized Differentiability

This section is devoted to the discussion of several concepts on generalized differentiability. We introduce the Bouligand and the slant derivative of a nonsmooth function, and we discuss the relationship between these two concepts. Further, we present a mean value theorem for Bouligand differentiable functions which also holds for semismooth functions.

Definition 2.1**.**

Let $X$ and $Y$ be two normed spaces, $D$ be a nonempty open set in $X$ and $J:D\subset X\rightarrow Y$ be a given mapping. For $x\in D$ and $h\in X$, if the limit

[TABLE]

exists, the function is said to be directionally differentiable. Further, $J^{\prime}(x)(h)$ is the directional derivative of $J$ at $x$ in the direction $h$.

Remark 2.2**.**

From here on, by making a small notation abuse, we denote by $F^{\prime}(u)$ the Fréchet derivative of $F$ at $u$, and by $F^{\prime}(u)(v)$ the directional derivative of $F$ at $u$ in the direction $v$.

Next, we define the concept of Bouligand differentiability and its relation to the semismoothnes concept. For further details, we refer the reader to [32, Ch. 2, Sec. 2.1].

Definition 2.3**.**

Let $D\subset\mathbb{R}^{n}$ be open and $J:D\rightarrow\mathbb{R}^{m}$ be Lipschitz continuos near $x\in D$, i.e., locally Lipschitz continuous. The set

[TABLE]

is called Bouligand-subdifferential (or B-subdifferential) of $J$ at $x$. Here, $D_{J}\subset D$ is the set of all $x\in D$ at which $J$ admits a Fréchet derivative $J^{\prime}(x)\in\mathbb{R}^{m\times n}$. The convex hull of the Bouligand-subdifferential $\partial J(x)=co(\partial_{B}J(x))$ is the Clarke’s generalized Jacobian of $J$ at $x$.

Definition 2.4**.**

Let $J:D\rightarrow\mathbb{R}^{m}$ be defined on the open set $D\subset\mathbb{R}^{n}$. $J$ is called B-differentiable at $x\in D$ if $J$ is directionally differentiable at $x$ and

[TABLE]

Proposition 2.5**.**

Let $D\subset\mathbb{R}^{n}$ be open and $J:D\rightarrow\mathbb{R}^{m}$ be a locally Lipschitz continuous function which is directionally differentiable at $x_{0}\in D$. Then, the function $J$ is B-differentiable at $x_{0}$.

Proof.

[29, Th.3.1.2 ] ∎

Let us notice that the B-derivative of a locally Lipschitz continuous function is its directional derivative.

Proposition 2.6**.**

Let $J:D\rightarrow\mathbb{R}^{m}$ be defined on the open set $D\subset\mathbb{R}^{n}$ and let $J$ be B-differentiable at x. Then, $J^{\prime}(x)(\cdot)$ is Lipschitz continuous, and, for every $h\in\mathbb{R}^{n}$, there exists $M\in\partial J(x)$ such that

[TABLE]

Proof.

[21, Sec. 8.2.1] ∎

Proposition 2.7**.**

Let $J:D\rightarrow\mathbb{R}^{m}$ be defined on the open set $D\subset\mathbb{R}^{n}$. Then, for $x\in D$ the following statements are equivalent:

$J$* is semismooth at $x$.* 2. 2.

$J$* is Lipschitz continuous near $x$, $J^{\prime}(x)(\cdot)$ exists, and*

[TABLE] 3. 3.

$J$* is Lipschitz continuous near $x$, $J^{\prime}(x)(\cdot)$ exists, and*

[TABLE]

Proof.

[32, Prop.2.7] ∎

Remark 2.8**.**

Let us notice that, from Proposition 2.5 we have that $J$ is B-differentiable if it is Lipschitz continuous near $x$ (locally Lispchitz continuous at $x$) and directionally differentiable at $x$. Then, from Proposition 2.7, we have that if $J$ is semismooth, $J$ is B-differentiable.

Proposition 2.9**.**

Let $J:D\rightarrow\mathbb{R}^{m}$ be defined on the open set $D\subset\mathbb{R}^{n}$ and let $J$ be semismooth in $x$. Then

[TABLE]

Proof.

[21, Th. 8.2] ∎

Next, we introduce the notion of slant differentiability [13, Sec.1] that will be used in this work.

Definition 2.10**.**

Let $X$ and $Y$ be Banach spaces, and $D$ be an open subset of $X$. A function $J:D\subset X\rightarrow Y$ is said to be slantly differentiable in $D$, if there exists a family of mappings $J^{\circ}:D\rightarrow\mathcal{L}(X,Y)$ such that

[TABLE]

for every $x\in D$. The function $J^{\circ}$ is called a slanting function for $J$ in $D$.

The previous definition was introduced in [13] as an adaptation of the definition of slant differentiability in Banch spaces stated in [5], where the family of linear operators $\{J^{\circ}(x+h)\}$ is required to be uniformly bounded in the operator norm. Also, in [13, pp. 868] the authors state that the notion of slant differentiability is a more general concept than the definition of semismoothness. In fact, the slanting functions are not required to be elements of the Clarke’s generalized Jacobian $\partial J(x+h)$. However, a single-valued selection $M(x+h)\in\partial J(x+h)$, with $x\in D$, is a slanting function $J^{\circ}$- in the sense of Definition 2.10 - if Proposition 2.7 (item 3) holds for $x\in D$ (see [13, Sec.1] for further details).

Proposition 2.11**.**

If $J$ is continuously differentiable in a neighbourhood of $x$, then

[TABLE]

Proof.

[32, Prop.2.2] ∎

We now present an important example of a semismooth function that will be useful in the subsequent sections.

Example 2.12**.**

Let $g>0$ be a constant. The mapping

[TABLE]

from $\mathbb{R}^{n}$ to $\mathbb{R}$ is semismooth on $\mathbb{R}^{n}$. Further, the slant derivative of this function is the characteristic function $\chi_{A_{\gamma}}(\vec{z})$ defined by

[TABLE]

where $A_{\gamma}:=\left\{\vec{z}:\gamma|\vec{z}|\geq g\right\}$

Proof.

[13, Sec.3, Lemma 3.1] ∎

Finally, we introduce the mean value theorem for B-differentiable functions.

Theorem 2.13**.**

(Mean value theorem for B-differentiable funcions.) Let $D\subset\mathbb{R}^{n}$ be an open convex set, $J:D\rightarrow\mathbb{R}^{m}$ be a B-differentiable function, and $x_{0},x_{1}\in D$. The function $\varphi:[0,1]\rightarrow\mathbb{R}^{m}$ defined by $\varphi(t)=J^{\prime}(x_{0}+t(x_{1}-x_{0}))(x_{1}-x_{0})$ is Lebesgue integrable and

[TABLE]

Proof.

See [29, Prop. 3.1.1]. ∎

3 Regularization and Discretization

The minimization problem (1) involves a convex non-smooth functional. In fact, the norm $|\nabla y|$ in the second term implies that the functional in (1) is not differentiable. In this context, we propose a regularization approach, based on a local Huber regularization procedure. Huber regularization has been used in previous contributions to approximate numerically several free boundary and nonsmooth problems with similar structure (see [12] and the references therein).

Let $\gamma>0$. We introduce the function $\psi_{\gamma}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ as follows:

[TABLE]

The function $\psi_{\gamma}$ corresponds to a local regularization of the Euclidean norm. Thanks to this procedure we obtain the following regularized optimization problem

[TABLE]

Theorem 3.1**.**

Let $1<p<\infty$ and $\gamma>0$. Then, problem (2) has a unique solution $u_{\gamma}\in W_{0}^{1,p}(\Omega)$. Also, the sequence $\{u_{\gamma}\}\subset W_{0}^{1,p}(\Omega)$ converges strongly in $W_{0}^{1,p}(\Omega)$ to the solution of problem (1), as $\gamma\rightarrow\infty$.

Proof.

See [12, Sec. 2]. ∎

3.1 Finite element approximation

Let us introduce the finite element approximation of problem (2). Let $\Omega_{h}$ be a given triangulation of the domain $\Omega$ , $n_{e}\in\mathbb{N}$ be the number of triangles $T_{i}$ such that $\bar{\Omega}_{h}=\displaystyle\cup_{i=1}^{n_{e}}T_{i}$ and $N$ be the number of nodes of the triangulation $\Omega_{h}$. For any two triangles, their closures are either disjoint or have a common vertex or a common edge. Finally, let $\{P_{j}\}_{j=1,\cdots,N}$ be the vertices (nodes) associated with $\Omega_{h}$. Taking this into account, we define

[TABLE]

where $\mathbb{P}_{1}$ is the space of continuous piecewise linear functions defined on $\Omega_{h}$. Then the following space

[TABLE]

is the finite-dimensional space associated with the triangulation $\Omega_{h}$.

Considering the previous analysis, the finite element approximation of (2) is formulated as follows

[TABLE]

Remark 3.2**.**

The finite element approximation of problems like (2) is restricted by the limited higher order regularity for the solution of the $p$-Laplacian (see [14]). Due to this fact, this kind of problems are usually approximated by continuous piecewise linear elements, which we also implement in this paper. In [1, 2, 24] the authors discuss optimal error estimates for sufficiently regular solutions, which can be obtained for specific data $f$ and $\Omega$. These results, however, are only valid for the $p$-Laplacian problem, i.e., when $g=0$. Since the solutions for variational inequalities of the second kind usually exhibit low global regularity, a deeper analysis is needed in order to obtain optimal error estimates for VIs. This, we consider, is beyond the scope of this paper.

The convergence analysis of the multigrid algorithm is based on the differentiability properties of the functional $J_{\gamma,h}(u)$. In order to analyze the regularity of this discrete functional, we decompose it as follows

[TABLE]

where

[TABLE]

$\mathcal{F}_{h}(u)$ is a functional associated with the discretized homogeneous Dirichlet problem for the $p$-Laplace operator, and it is known to be a twice Fréchet-differentiable and strictly convex functional with Fréchet derivative $\mathcal{F}^{\prime}_{h}(u_{h})$ given by

[TABLE]

and

[TABLE]

See ([1, 2, 9, 12]) for further details.

Proposition 3.3**.**

Let $1<p<\infty$. The functional $J_{\gamma,h}(u_{h})$ is differentiable with

[TABLE]

Furthermore $J^{\prime}_{\gamma,h}(u_{h})$ is semismooth and its slant derivative is given by

[TABLE]

Proof.

Let us start by analyzing the functional $\mathcal{G}_{\gamma,h}(\nabla u_{h})$. It is known that $\mathcal{G}_{\gamma,h}$ is differentiable (see [12, Sec. 2.2]), and moreover, we know that

[TABLE]

where

[TABLE]

By using the $\max$ function, we can rewrite $\mathcal{G}^{\prime}_{\gamma,h}(\nabla u_{h})(v_{h})$ in the following way.

[TABLE]

Next, from (5), it follows that

[TABLE]

which, thanks to (6) and (10), implies that

[TABLE]

The second derivative of $J_{\gamma,h}(u_{h})$ does not exist. In fact, the functional $\mathcal{G}^{\prime}_{\gamma,h}(\nabla u_{h})$ is not differentiable since this functional involves the $\max$ function. Fortunately, the $\max$ function is slantly differentiable when defined in finite dimensional spaces (see Example 2.12). Thus, we can calculate the slant derivative of $\mathcal{G}^{\prime}_{\gamma,h}(\nabla u_{h})$, denoted by $(\mathcal{G}^{\prime}_{\gamma,h})^{\circ}(\nabla u_{h})$, as follows.

$|\nabla u|\geq\frac{g}{\gamma}$: Here, we have that

[TABLE]

where $\chi_{A_{\gamma,h}}$ is the slant derivative of the function $\max(g,\gamma|\nabla u_{h}|)$.

$|\nabla u|<\frac{g}{\gamma}$: Here, we have that

[TABLE]

where $I_{\gamma,h}:=\Omega_{h}\setminus A_{\gamma,h}$.

Then, the slant derivative of $\mathcal{G}^{\prime}_{\gamma,h}(\nabla u_{h})$ reads as follows

[TABLE]

On the other hand, from (7) we have that

[TABLE]

Hence, from (11) and Proposition 2.11, we can state that

[TABLE]

which, thanks to (12), (13) and (14), yields that

[TABLE]

∎

Proposition 3.4**.**

The slanting function $(J^{\prime}_{\gamma,h})^{\circ}(u_{h})$ is positive definite.

Proof.

Following the decompostition presented in (14), we know that

[TABLE]

It is well known that $\mathcal{F}$ is a strictly convex functional, which implies that

[TABLE]

Next, let us recall the expression $(\mathcal{G}^{\prime}_{\gamma,h})^{\circ}(u_{h})(w_{h},w_{h})$ given by

[TABLE]

Applying Cauchy-Schwarz to the right hand side in (18), we have that

[TABLE]

which implies that

[TABLE]

Then, $(J^{\prime}_{\gamma,h})^{\circ}(u_{h})(w_{h},w_{h})>0$. ∎

4 The Multigrid for Optimization (MG/OPT) Method

In this section, we present the multigrid for optimization (MG/OPT) algorithm for solving the regularized and discretized optimization problem (4). The MG/OPT method was introduced as an efficient tool for large scale optimization problems (see [25, 23]). The idea of the algorithm is to take advantage of the solutions of problems discretized in coarse meshes to optimize problems in finer meshes. The efficient resolution of coarse problems provides a way to calculate search directions for problems discretized in finer levels.

In order to present the algorithm, we shall introduce the following preliminaries. Let $\Omega$ be a given bounded domain. A standard procedure to generate a sequence of triangulations on $\Omega$ is to define a coarse mesh and then refine it several times until getting the sequence $\{\Omega_{k}\}_{k=0,\ldots,m}$, i.e., each $\Omega_{k}$ is obtained from $\Omega_{k-1}$ by regular subdivision. This procedure joins the edge midpoints of any triangle in mesh $\Omega_{k-1}$ by edges, and forms the new triangles of $\Omega_{k}$ (see Figure 1). Each node in $\Omega_{k-1}$ is also a node in $\Omega_{k}$, and every node belonging to $\Omega_{k}$ but not to $\Omega_{k-1}$ is the midpoint of an edge in $\Omega_{k-1}$. Let us denote by $N^{(k)}$ the number of nodes associated to each $\Omega_{k}$. Thus, the nodes of $\Omega_{k-1}$ are the first $N^{(k-1)}$ nodes of $\Omega_{k}$ (see Figure 2).

When working with domains with curved boundaries, the refinement procedure described above is slightly different in the boundary edges. This is the case, for instance, of the domain depicted in Figure 2. During the refinement, instead of generating a new midpoint in the nearest triangulation edges to the boundary, each edge is replaced by two edges that intersect at the midpoint of the curved-segment of the boundary, i.e., the boundary is bisected by the two new edges added. Consequently, the resulting triangulation only covers the domain $\Omega$ approximately, introducing a source of error [11, Ch. 4, pp. 93]. Since a non polygonal domain can not be triangulated exactly, approximating a domain with a curved boundary is a matter under research. One interesting technique to analyze this problem is the isoparametric method which uses finite elements with curved edges [11, Ch.4, Sec. 4.7]. This approach seems to be the most appropriate for problems of fluid mechanics, which is why it will be considered in future contributions.

The multigrid approach involves two transfer operators. As we are working with a set of meshes and the algorithm runs at each level of discretization, we need to transfer information among the different grids. Then, we intoduce the fine-to-coarse grid transfer operator, $I_{k}^{k-1}$, and the coarse-to-fine grid transfer operator, $I_{k-1}^{k}$.

The latter operator transfers information from a coarse mesh $\Omega_{k-1}$ to a finer mesh $\Omega_{k}$. It is also called the prolongation operator. On the other hand, $I_{k}^{k-1}$ or restriction operator, transfers information from a fine grid to a coarse one.

In this paper we use the mesh data structure and the operators implemented by M. S. Gockenbach in order to obtain $I_{k}^{k-1}$ and $I_{k-1}^{k}$ (see [11, Ch. 6, Ch. 13]). In what follows, the ideas stated in the previous reference are outlined for the prolongation and restriction procedures in a coarse grid $\Omega_{k-1}$ and a finer grid $\Omega_{k}$. For instance, in Figure 2 we have the first stage of refinement in a disk domain with grids $\Omega_{0}$ and $\Omega_{1}$ .

In order to transfer information from the nodes in a coarse mesh $\Omega_{k-1}$ (e.g. $\Omega_{0}$ in Figure 2) to a refined mesh $\Omega_{k}$ (e.g. $\Omega_{1}$), we first copy the nodal values of the nodes in $\Omega_{k-1}$ to the corresponding nodes in $\Omega_{k}$ (in our example they correspond to the nodes with indexes $1,2,3,4$ and $5$). By construction, we know that these are the first $N^{(k-1)}$ ($N^{(0)}=5$) nodes of $\Omega_{k}$, and,the midpoints complete the $N^{(k)}$ ($N^{(1)}=13$) nodes of the mesh $\Omega_{k}$ (in $\Omega_{1}$ they correspond to the nodes with indexes $6,7,\cdots,13$). Then, since the values for nodes with index $j=1,\cdots N^{(k-1)}$ are already given, we need to compute the nodal values for the midpoints, i.e., for nodes with index $j=N^{(k-1)}+1,\cdots,N^{(k)}$. Let $V(\textbf{P})_{k}=\{V(P_{j})\}_{j=1,\cdots,N^{(k)}}$ be the vector containing the nodal values associated to $\Omega_{k}$ The prolongation operator computes the values to the midpoints as follows:

[TABLE]

where the nodes $P_{end(j,1)}$ and $P_{end(j,2)}$ are the endpoints of the edge in $\Omega_{k-1}$ for which $P_{j}$ is the midpoint. These endpoints are called node parents and, by construction, are unique for each node.

Now, to transfer information from the nodes in $\Omega_{k}$ to $\Omega_{k-1}$ (e.g. $\Omega_{1}$ to $\Omega_{0}$ in Figure 2) the restriction operator first copy the values of the nodes in $\Omega_{k}$ to the corresponding nodes that are also in $\Omega_{k-1}$. Then, for $j>N^{k-1}$, the operator perform the following computation

[TABLE]

The mesh data structure implemented allows us to copy and extract nodes from any specific discretization. We refer the reader to see [11, Ch. 13, sec. 13.2.1] for a detailed explanation of the transfer operators.

In multigrid schemes is standard to assume that

[TABLE]

where $c$ is a constant. In our case, if we compute the products $I_{k-1}^{k}V(\textbf{P})_{k-1}=V(\textbf{P})_{k}$ and $I_{k}^{k-1}V(\textbf{P})_{k}=V(\textbf{P})_{k-1}$, the restriction and prolongation operators satisfy the condition (see [11, Sec. 13.2.1, p. 294])

[TABLE]

Now that we have introduced the interpolation operators, we are ready to discuss the MG/OPT algorithm for problem (4). The MG/OPT method corresponds to a nonlinear programming adaptation of the full approximation storage (FAS) scheme (see [4, 30]). The multigrid subproblems arising from the different discretization levels are nonlinear optimization problems [23]. Then, MG/OPT is related to different optimization techniques ranging from the gradient method to quasi Newton methods to solve the problems at each level. The multigrid for optimization approach requires mild conditions regarding the underlying optimization algorithm. Mainly, this algorithm needs to be globally convergent. However, it is important to highlight that at each level of discretization we need to find a sufficiently accurate solution for the minimization subproblem. Then, the selection of the underlying optimization algorithm is not trivial and depends on the inner characteristics of the optimization problem.

In the present problem, since the $p$-Laplacian is involved in the objective functional $J_{\gamma}$, we have to consider that its finite element approximation (4) results in a nonlinear and possibly degenerate finite dimensional problem (this may be the case when $p<2$, see [14] and the references therein). Also, the functional $J_{\gamma,h}$ involves a semismooth function. These facts need to be taken into account when proposing the underlying optimization algorithms. In our case, we propose a class of preconditioned descent algorithms, designed specifically for $p$-Laplacian type problems. For the convenience of the reader, in the next section we describe briefly this algorithm as is implemented in our computational results.

As we mentioned before, the main idea of the MG/OPT algorithm is to use coarse problems to generate, recursively, search directions for finer problems. Then, a line search procedure, along with the underlying optimization algorithm is used to improve the solution at each level of discretization.

In what follows we present the MG/OPT algorithm. The underlying optimization algorithm will be denoted by $S_{opt}$ inside the multigrid approach. The initial discretized problem is given on the finest grid. To facilitate the implementation of the algorithm, the MG/OPT scheme is presented in a recursive formulation. Hence, we introduce the following slightly different notation for the optimization problem

[TABLE]

We set $\hat{f}_{k}=0$ at the finest level $k=m$. Then, $\hat{J}_{\gamma,k}$ corresponds to the functional $J_{\gamma,h}$ introduced in problem (4), discretized at each level $k=0,\ldots,m$. Therefore, we replace the subscript $h$ by $k$. Hereafter, the same treatment is given to $u_{k}$, which stands for $u_{h}$, discretized at each level $k$.

Remark 4.1**.**

Since $J_{\gamma,k}$ is Fréchet differentiable and defined in finite dimensional spaces, we can associate the derivative $J^{\prime}_{\gamma,k}(u_{k})$ with the gradient $\nabla J_{\gamma,k}(u_{h})$ as follows (see [10, Sec.1.1, p. 8])

[TABLE]

In what follows, we will use this representation.

Summarizing and taking into account Remark 4.1, the algorithm reads as follows.

The algorithm presented above contemplates one iteration of a V-cycle initialized with a rough estimate of the solution on the finest grid.

5 Convergence Analysis

In this section, we discuss the convergence properties of Algorithm 1. Following [23, 25], we can state that the global convergence of the underlying optimization algorithm ensures global convergence of the MG/OPT method. This comes from the fact that if we have an approximate solution (given by the underlying optimization algorithm) at each discretization level, the algorithm generates search directions for problems discretized on finer meshes. Once we have the descent direction, a line search procedure is used to improve the solution at each finer problem.

In the classical convergence analysis of the MG/OPT methods, the three following conditions are critical (see [23, 25, 31]).

The discretized objective functional is convex at each level of discretization. 2. 2.

The subproblems

[TABLE]

are solved accurately enough. 3. 3.

The transfer operators satisfy the standard condition $I_{k}^{k-1}=c\left(I_{k-1}^{k}\right)^{T}$.

These conditions are helpful to prove that the search direction provided by the MG/OPT algorithm is indeed a descent direction. For instance, the convexity condition is key to prove that the Hessian is positive definite. However, in our case, the classical Hessian does not exist. Thus, we use the slant differentiability of the functional $J_{\gamma,h}$ to obtain positive definiteness, see Proposition 3.4.

Next, we comment on the assumptions on our problem. We are considering the subproblems (21) instead of the subproblems $\displaystyle\min_{u_{k-1}}\left(\hat{J}_{\gamma,k-1}(u_{k-1})-\hat{f}^{\top}u_{k-1}\right)$ presented in Algorithm 1. The artificial term $I_{k}^{k-1}\hat{f}_{k}$ is simply introduced in order to facilitate the recursive implementation of the method. In fact, at the very beginning of the algorithm, on the finest mesh, this term is set to zero. In the succeeding iterations, $I_{k}^{k-1}\hat{f}_{k}$ corresponds to the recursive sum of the previous restricted $\tau_{k-1}$. Thus, in order to make the subsequent analysis easier, we analyze the subproblems (21). Also, since problem (4) is strictly convex in $V_{h}^{0}$, without loss of generality, subproblems (21) are strictly convex in $V_{k}^{0}$, for $k=0,\ldots,m$.

Since we perform a few iterations of a suitable globally convergent optimization algorithm ($S_{opt}$), we ensure that the subproblems $\min_{u_{k-1}}\left(J_{\gamma,k-1}(u_{k-1})-\tau_{k-1}^{\top}u_{k-1}\right)$ are solved accurately enough. Finally, (20) yields that $I_{k}^{k-1}=\left(I_{k-1}^{k}\right)^{T}$.

Let us recall that the search direction for the MG/OPT algorithm is denoted by $e$ and search directions of the underlying optimization algorithm, (inside the MG/OPT loop) are denoted by $w_{k}$.

Proposition 5.1**.**

The search direction $e=I_{k-1}^{k}(\widetilde{u}_{k-1}-u_{k-1}^{\nu_{1}})$ is a descent direction for all $k=1,\ldots,m.$

Proof.

We have to prove that

[TABLE]

From this point, for the readability of the proof, we drop the subscript $\gamma$. First note that, if we solve

[TABLE]

exactly, then

[TABLE]

Since we are solving the problem inexactly (but accurately enough), we then have that,

[TABLE]

for some $z$ as small as the algorithm accuracy allows for. From Algorithm 1 we have that

[TABLE]

Hence, we can rewrite (23) as follows

[TABLE]

This expression implies that

[TABLE]

where

[TABLE]

Next, let us focus on the first term in the right-hand side of (25). Hereafter, we use the notation $\nabla J_{k-1}=H_{k-1}$. Then, we have that

[TABLE]

We know, from Proposition 3.3, that $H_{k-1}$ is a semismooth function, which, thanks to Remark 2.8, implies that $H_{k-1}$ is B-differentiable. Thus, Theorem 2.13 yields that

[TABLE]

where $H_{k-1}^{\prime}(u_{k-1}^{\nu_{1}}+t(w_{k-1}))(w_{k-1})$ stands for the directional derivative of the operator $H_{k-1}$ at $(u_{k-1}^{\nu_{1}}+t(w_{k-1}))$ in the direction $w_{k-1}$. Furthermore, $H_{k-1}^{\circ}(\widetilde{u}_{k-1})(w_{k-1},w_{k-1})$ is given by (9). Following (14) and Proposition 3.4, we know that $H_{k-1}^{\circ}(\widetilde{u}_{k-1})$ is definite positive. Consequently, there exists a constant $c$ such that

[TABLE]

Next, let us focus on the last term on the right hand side of (27). By using the Cauchy-Schwarz inequality and the fact that $H_{k-1}^{\prime}(u_{k-1}^{\nu_{1}}+t(w_{k-1}))(w_{k-1})$ is Lebesgue integrable, we obtain that

[TABLE]

Moreover, since $H_{k-1}$ is semismooth, from Propositions 2.7 and 2.9 we have that

[TABLE]

Hence, for an arbitrary $\epsilon>0$, it holds that

[TABLE]

Thus, from (29) we have that

[TABLE]

Finally, by taking $\epsilon<\frac{c}{2}$ and considering (27), (28) and (30), we conclude that

[TABLE]

Note that $w_{k-1}=\widetilde{u}_{k-1}-u^{\nu_{1}}_{k-1}\neq 0$. Consequently, we have that

[TABLE]

In order to prove that $e$ is a descent direction, we still need to prove that the third term of the right hand side in (25) satifies that

[TABLE]

Note that $\widetilde{u}_{k-1}$ is the solution of the problem

[TABLE]

Therefore,

[TABLE]

which is equivalent to

[TABLE]

since the optimization algorithm was initialized with $u_{k-1}^{\nu_{1}}$. Now, using the mean value theorem for differentiable functionals we have that

[TABLE]

for some $\xi$ between $\widetilde{u}_{k-1}$ and $u_{k-1}^{\nu_{1}}$. Hence, from the inequality (33) and equation (34) we have that

[TABLE]

which implies that

[TABLE]

Next, by approximating $\xi\approx\widetilde{u}_{k-1}$, from (23) and (35), we obtain that

[TABLE]

Summarizing, (25), (32) and (36) imply that

[TABLE]

and we can conclude that $e$ is a descent direction. ∎

Finally, thanks to the previous results, we can state and prove the following theorem of convergence for the MG/OPT Algorithm 1.

Theorem 5.2**.**

Suppose that the following hypotheses are satisfied:

• •

The optimization algorithm, $S_{opt}$, applied to an optimization problem of any resolution, is globally convergent, i.e.,

[TABLE]

• •

At least one of the parameters $\nu_{1}$ or $\nu_{2}$ is positive.

• •

The search direction $e=I_{k-1}^{k}(\widetilde{u}_{k-1}-u_{k-1}^{\nu_{1}})$ is a descent direction.

Then the MG/OPT algorithm is globally convergent in the sense of (37).

Proof.

Let us start by noticing that if $\nu_{1}$ or $\nu_{2}$ are positive, at least one iteration of the optimization algorithm is performed at every cycle of MG/OPT. Thus, an approximate solution at each level of discretization $k$ is obtained. Since the search direction $e=I_{k-1}^{k}(\widetilde{u}_{k-1}-u_{k-1}^{\nu_{1}})$ is a descent direction, the approximate solution given at each level improves at every cycle of MG/OPT, i.e., the functional value $\hat{J}_{\gamma,k}(u_{k})$ decreases at each cycle after the solution update $u_{k}^{\nu_{1}+1}=u_{k}^{\nu_{1}}+\alpha_{k}e$. Consequently, as the underlying optimization algorithm is globally convergent, the multigrid optimization algorithm is globally convergent. ∎

Remark 5.3**.**

The convergence of MG/OPT methods depends on the behavior of the underlying algorithms. In fact, usually the performance of the selected optimization method is the same that is verified for the MG/OPT algorithm (see [3, 23, 25, 31]). Moreover, it is common that for optimization methods designed for nonsmooth problems, only global convergence results can be established. In our particular selection for the underlying algorithm, the global convergence of the preconditioned descent method is guaranteed [12]. However, the rates of this convergence are not investigated. To the best of our knowledge, the use of semismooth derivatives, in the context of multigrid methods, is a novel perspective. In the present case, we need to analyse the behavior of the “slant Hessian” $(J^{\prime}_{\gamma,h})^{\circ}$ in the context of the underlying algorithms. This idea is currently under study and it will be developed in future contributions.

6 Implementation

6.1 Optimization algorithm

In this section, we briefly discuss the preconditioned descent algorithm proposed in [12], which is proposed for the underlying optimization algorithm within the MG/OPT algorithm.

Generally speaking, a descent method starts with an initial point $u_{0}$ and, with information of first order, the algorithm finds directions that lead us to the minimum of the objective functional. Also, the method must find the length of the step, $\alpha_{r}$, along the chosen direction, $w_{r}$. The basic idea consists in finding $\alpha_{r}$ and $w_{r}$ such that:

[TABLE]

in every iteration of the method.

In the gradient method, for instance, the search direction $w_{r}$ is determined by

[TABLE]

On the other hand, for the preconditioned descent algorithm proposed in [12], the search direction $w_{r}$ is determined by solving the following variational equation

[TABLE]

where the form $P_{r}:V_{h}^{0}\times V_{h}^{0}\rightarrow\mathbb{R}$ is chosen as a variational approximation of the $p$-Laplacian.

The algorithm reads as follows

The global convergence and the stability of this Algorithm is guaranteed both in finite and infinite dimension settings. For a deeper analysis of this algorithm, we refer the reader to [12].

6.2 Line search technique

In this section, we describe the line search algorithm which will be used in the implementation of Algorithm 1. This algorithm uses polynomial models of the objective functional for backtracking, and it was originally proposed in [7, Sec. 6.3.2].

The algorithm reads as follows

If we set

[TABLE]

then the quadratic model $m_{2}$ is given by

[TABLE]

while the cubic model $m_{3}$ is given by

[TABLE]

where

[TABLE]

and $\alpha_{p}$ and $\alpha_{2p}$ are the last two previous values of $\alpha_{k}$. For further details and examples see [7, pp. 126-129] and the references therein.

7 Applications to the Numerical Simulation of Viscoplastic Flow

In this section, we discuss the application of the MG/OPT algorithm to the numerical simulation of the steady flow of viscoplastic fluids. These materials are characterized by the existence of a yield stress [6, 8, 16, 15]. This implies that the viscoplastic material exhibits no deformation if the shear stress imposed does not exceed the yield stress, i.e., it behaves like an ideal rigid solid. If the shear stress overpasses the yield stress, the material will deform like a nonlinear viscous fluid in most of the cases. In fact, Herschel - Bulkley and Casson fluids present a nonlinear stress-shear rate, while Bingham fluids behave like a viscous fluid with linear stress-shear rate (see Figure 3). Summarizing, the existence of the yield stress causes rigid zones and yielded zones in the flow.

In this work we focus on the stationary and laminar flow of a viscoplastic fluid in a cylindrical pipe under the effect of a drop in pressure. We consider three classic models for these fluids: Herschel-Bulkley, Bingham and Casson. This kind of flow is a simplified problem in which we can assume that all velocity fields only have a non-zero component in the axial direction. Therefore, by assuming that the velocity fields vanish on the boundary of the pipe (the so called adhesion condition) and that $f$ represents the constant pressure drop, it is well known that the velocity field across the cross-section of the pipe solves the following optimization problem.

[TABLE]

where

[TABLE]

This variational formulation is motivated by the necessity of representing the free surface that separates the regions in which the material has yielded from those in which it behaves like a rigid solid. Regarding more general problems, like the flow in 2D and 3D geometries, the variational formulation also leads to optimization problems formulated in divergence-free spaces. This fact suggests that the present methodology can be generalized to these more challenging problems. We refer the reader to [8, 12, 16] for a more detailed explanation of the variational approach to the flow problems.

In order to validate our results, we introduce the theoretical velocity distribution, in cylindrical coordinates, for a circular pipe flow. Here $r=r_{0}=2g$ (See [16]).

• •

Herschel-Bulkley

[TABLE]

with $\beta=\frac{1}{p-1}$.

• •

Bingham

[TABLE]

• •

Casson

[TABLE]

Hereafter, we discuss the performance of the MG/OPT algorithm when applied to the numerical solution of (38). All the numerical experiments in this paper are implemented in MATLAB (R2015a) and run on an Intel Core i5 processor with 2.5 GHz. The MG/OPT algorithm was implemented in a V-cycle scheme with $\nu_{1}$ and $\nu_{2}$ as the pre and post optimization iterations. In all experiments the right hand side was set $f=1$. We initialize the algorithms with the solution of the Poisson problem $-\Delta u_{h}=f$. We performed the numerical experiments in the unit circle domain. Also, we compare the performance of the MG/OPT algorithm with the performance of the preconditioned descent algorithm when solving the same problem in the finest grid.

Regarding the stopping criteria, we stop the algorithms 1 and 2 as soon as the expression $\|\nabla J_{k}(u_{k})\|$ is reduced by a factor of $10^{-7}$.

7.1 Herschel-Bulkley

The Herschel-Bulkley model can be seen as a power-law model with plasticity [12, 16]. The introduction of the parameter $p$, known as the flow index, allows the model to represent the behaviour of several viscoplastic fluids. In fact, if $1<p<2$, the material exhibits a pseudoplastic or shear-thinning behaviour. On the other hand, if $p>2$ the fluid behaves as a shear-thickening material (see Figure 3). This versatility makes the Herschel-Bulkley model to be widely used to simulate several materials including liquid foams, whipped cream, fluid foods, silly putty and some polymers [6].

Experiment 1 (case $1<p<2$)

In this experiment we compute the solution for problem (38) with $p=1.75$ and different values of $g$, ranging from $g=0$ to $g=0.4$. We fix the preconditioned descent algorithm (see [12]) as the underlying optimization algorithm $S_{opt}$, with set $\epsilon=10^{-6}$. MG/OPT V-cycles are carried out in 5 grids, with $8321$ nodes in the finest grid and $41$ nodes in the coarsest one.

In Table 1, we compare the MG/OPT algorithm and the preconditioned descent algorithm. Following [25], we stop the MG/OPT Algorithm 1 at a given number of iterations and compare the numerical performance with the preconditioned descent algorithm. In particular, we compare the error $Err_{s}$, which is defined as the difference between the current solution and a highly accurate solution, measured in the discrete $L_{2}$-norm. The accurate solution is obtained by running the globally convergent preconditioned descent algorithm on the finest mesh. We report $Err_{s}$ in the third and seventh columns for MG/OPT and the preconditioned algorithm, respectively. Also, since we know the analytical solution for the velocity distribution along the radio in the circular pipe (39), in the fifth and eight columns of Table 1 we present the error for the constant plug flow velocity , i.e., $Err_{pf}=|u_{\text{a}}-u_{\text{num}}|$, where $u_{\text{a}}$ is the analytical plug flow velocity calculated with (39) and $u_{\text{num}}$ stands for the plug flow velocity calculated numerically with both algorithms. With $g=0$, MG/OPT performs less iterations than the preconditioned descent algorithm. However, since at each iteration or V-cycle of MG/OPT the algorithm performs $\nu_{1}+\nu_{2}=4$ optimization steps at the finest grid, the total number of optimization steps in the finest grid is not reduced substancially in comparison with the preconditioned descent algorithm. In the subsequent cases with a larger $g$, MG/OPT iterations are about one seventh of the preconditioned descent algorithm iterations. Also, in comparison with the underlying optimization algorithm, MG/OPT reduces by almost half the number of optimization steps in the finest grid at every stage of the computation, with $g=0.2$. This behaviour is replicated in the advanced stages of MG/OPT with $g=0.4$.

The resulting velocity fields $u_{h}$ are displayed in Figure 4. The regions with constant velocities in the center of the pipe corresponds to the plug flow velocity, where the material behaves like a rigid body. These regions, of course, do not appear with $g=0$.

7.2 Bingham

Bingham fluids are viscoplastic materials that can be seen as a particular case of the Herschel-Bulkley model when $p=2$. The main characteristic of Bingham fluids is that when the shear stress exceeds the yield stress, the material exhibits a close to linear stress-strain relation (see Figure 3).

Experiment 2

The following experiment was carried out in a disk domain with the same characteristics as the previous experiments. We consider the following parameters: $g=0.4$ and $\gamma=10^{3}$. In this experiment, we execute a comparison in the same terms as in the previous experiment. In Table 2, MG/OPT with $\nu_{1}=\nu_{2}=2$ does not reduce the optimization steps in the finest grid (each iteration of MG/OPT performs $4$ optimization steps in the finest grid) and $Err_{s}$ hardly decreases. Due to this unexpected behaviour, the number of optimization steps in the post-optimization procedure was increased and decreased in the pre-optimization procedure by setting $\nu_{1}=1$ and $\nu_{2}=3$. This new computation is tabulated In Table 3. In this case, the optimization steps in the finest grid are reduced substantially. The numerical experience in this experiment tells us that the restriction procedure introduce errors to the solution at each level. Then, these errors are reduced more efficiently by increasing the number of optimization steps in the post optimization phase.

[TABLE]

The resulting velocity field is displayed in Figure 5. Since the yield stress $g=0.4$ is high, the plug flow covers a large part of the cross section of the pipe.

7.3 Herschel-Bulkley: case $p>2$

We now analyze the behaviour of the MG/OPT algorithm in the case $p>2$, which usually gives rise to instabilities in the performance of the numerical algorithms (see [1, 2, 9, 14]).

Experiment 3

In the following experiment we present the performance of the MG/OPT algorithm for Herschel-Bulkley with $p=5$ and $g=0.1$. When $p$ increases, the optimization problem (38) is difficult to solve [14]. Additionally, since we restrict the approximate solution to coarser grids, at the coarsest mesh the restriction procedure fails in approximating efficiently the plug flow. Hence, in this case we only consider MG/OPT in a V-cycle with 3 grids instead of 5 grids (in order to avoid several restriction procedures). Therefore, we work with $2113$ nodes in the finest grid and $145$ nodes in the coarsest one. The preconditioned descent algorithm [12] is fixed as the underlying optimization algorithm $S_{opt}$. With the previous setting, in Table 4 the MG/OPT algorithm (with $\nu_{1}=\nu_{2}=2$) does not reduce the optimization steps in the finest grid in the first two cycles. However, in the third cycle MG/OPT performs one quarter of the optimization steps of the preconditioned descent algorithm in the finest grid.

The resulting velocity field is displayed in Figure 6. Now we are in the case of a shear-thickening material. Since the shear stress transmitted by a fluid layer decreases toward the center of the pipe, the velocity takes a conical form.

7.4 Casson

The Casson model is a viscoplastic model that was first developed for modeling printing inks. It has also been used to model food flow behaviour such as chocolate and cocoa products [27], and has been applied to biorheology models like hemodynamics and viscometric flows [33]. The Casson model turns out to be special since it involves the sum of the terms $\frac{1}{2}\displaystyle\int_{\Omega_{h}}|\nabla u_{h}|^{2}\,dx+\frac{4}{3}\sqrt{g}\int_{\Omega_{h}}|\nabla u_{h}|^{\frac{3}{2}}\,dx$. The first term corresponds to the Bingham structure and the second one to the Herschel-Bulkley structure with $p=\frac{4}{3}<2$. Thus, we decide to use the preconditioned descent algorithm with $\epsilon=10^{-6}$.

Experiment 4

In the following experiment we present the performance of the MG/OPT algorithm for the Casson problem with parameter $g=0.2$. In this experiment we show the behaviour of MG/OPT with a full multigrid (FMG) scheme [4]. Here, our aim is to show the time reduction when solving the Casson problem with MG/OPT and without it. We compare this scheme with a different number of grids. In Table 5, the first row corresponds to the solution computed by the preconditioned descent algorithm in a circle domain in a fine grid with $8321$ nodes. We take as reference the number of iterations and the execution time of this algorithm i.e., the CPU time until the error $Err_{s}$ reaches a tolerance of $1e-07$. The subsequent rows correspond to the solution computed by the MG/OPT $(\nu_{1}=\nu_{2}=2)$ algorithm. For each row we present the number of grids used and the number of nodes in the coarsest grid. At each row, a new grid is added to the scheme.

Since we know the analytical solution (41) in the circular pipe, we present the error $Err_{pf}$ for the plug flow velocity in the fifth column. FMG performs $1$ V-cycle with $\nu_{1}=\nu_{2}=2$ optimization steps in the finest grid, totally we have 4 finest grid optimization steps whereas the preconditioned descent algorithm performed 10 iterations in order to achieve the same error. Let us notice that, as more grids are added to the scheme, we achieve higher CPU time savings in all cases. With 5 grids the time reduction achieved was around 40%.

The MATLAB codes are available in www.modemat.epn.edu.ec/$\sim$sgonzalez/Publications.html.

8 Conclusions and Outlook

We proposed and analyzed a multigrid for optimization (MG/OPT) algorithm for the numerical solution of a class of quasilinear variational inequalities of the second kind. We analyzed the variational inequality via the minimization of the associated energy functional. First, we regularized the non-differentiable part of the energy by using a Huber regularization approach. Next, we proposed a finite element discretization for the problem, and we analyzed the differentiability of the discretized functional. In particular, we proved that the Jacobian of the discretized functional is slantly differentiable. The MG/OPT algorithm was presented and all of the involved transfer operators analyzed. The convergence of the MG/OPT algorithm was established by using the mean value theorem for slantly differentiable functions and the global convergence of the underlying optimization algorithm. The main issues regarding the implementation of the algorithm were explained, and we described the type of global convergent deepest descent methods used as underlying optimization algorithms. We showed that several classical models for viscoplastic flow correspond to the class of variational inequalities under study. Therefore, we focussed the numerical experiments on this kind of problems. Particularly, we computed the solution for the Herschel-Bulkley, Bingham and Casson models. In all the experiments presented, we observed CPU-time savings, especially when working on the finest meshes. This showed that the MG/OPT algorithm is indeed an efficient tool for dealing with large scale problems.

In order to continue this research, we consider that the study of a more general class of variational inequalities is an interesting perspective. Since the functional $J_{\gamma,h}$ has a second slant derivative, the development of generalized Newton-type methods as smoothers for the FAS scheme looks like a promising field of research. Finally, the analysis of more challenging problems, such as the $p$-Stokes problem, problems coming from glaceology, geophysics and hemodynamics is of great interest.

Acknowledgements

We thank Prof. Juan Carlos De los Reyes (MODEMAT-EPN) for his keen observations in the theory of semismooth functions. We also thank the anonymous referees for many helpful comments which lead to a significant improvement of the article.