Gradient flows for the p -Laplacian arising from biological network models: A novel dynamical relaxation approach

Abstract

We investigate a scalar partial differential equation model for the formation of biological transportation networks. Starting from a discrete graph-based formulation on equilateral triangulations, we rigorously derive the corresponding continuum energy functional as the Γ-limit under network refinement and establish the existence of global minimizers. The model possesses a gradient-flow structure whose steady states coincide with solutions of the p-Laplacian equation. Building on this connection, we implement finite element discretizations and propose a novel dynamical relaxation scheme that achieves optimal convergence rates in manufactured tests and exhibits mesh-independent performance, with the number of time steps, nonlinear iterations, and linear solves remaining stable under uniform mesh refinement. Numerical experiments confirm both the ability of the scalar model to reproduce biologically relevant network patterns and its effectiveness as a computationally efficient relaxation strategy for solving p-Laplacian equations for large exponents p.

1. Introduction

In this paper, we study the scalar partial differential equation system modeling the formation of biological transportation networks,

$$\begin{array}{ccccc}\multicolumn{1}{c}{-\mathrm{\nabla }\cdot [(c+r)\mathrm{\nabla }u]}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{S,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(1)

$$\begin{array}{ccccc}\multicolumn{1}{c}{\frac{\mathrm{\partial c}}{\mathrm{\partial t}}-|\mathrm{\nabla }u{|}^2+\nu c^{\gamma -1}}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{0,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(2)

for the scalar pressure field $u=u(t,x)\in \mathrm{ℝ}$ of the fluid transported within the network and the scalar conductance (permeability) $c=c(t,x)\in \mathrm{ℝ}$ , with metabolic constant $\nu >0$ and metabolic exponent $\gamma >0$ . We assume the validity of Darcy’s law for slow flow in porous media, connecting the flux q of the fluid with the gradient of its pressure via the action of the permeability $\mathrm{\mathbb{P}}[c]$ ,

$$\begin{array}{ccc}\multicolumn{1}{c}{q=-\mathrm{\mathbb{P}}[c]\mathrm{\nabla }u.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

The total permeability is assumed to be of the form $\mathrm{\mathbb{P}}[c]:=c+r$ , where the scalar function $r(x)\ge 0$ describes the isotropic background permeability of the medium. The source term $S(x)$ in the mass conservation equation (1) is to be supplemented as an input datum and is assumed to be independent of time. We pose (1), (2) on a bounded domain $\mathrm{\Omega }\subset {\mathrm{ℝ}}^d$ , $d\ge 1$ , with smooth boundary $\partial \mathrm{\Omega }$ .

For (1), we prescribe the homogeneous Neumann boundary condition

$$\begin{array}{ccc}\multicolumn{1}{c}{\nu \cdot \mathrm{\nabla }u=0\mathrm{on}\partial \mathrm{\Omega },}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(3)

where $\nu =\nu (x)\in {\mathrm{ℝ}}^d$ is the outer normal vector to the boundary $\partial \mathrm{\Omega }$ . We note that (2) is a system of ordinary differential equations (ODEs) (coupled through the term $|\mathrm{\nabla }u{|}^2$ ) and, therefore, no boundary condition can be prescribed for $c=c(t,x)$ .

We prescribe the initial condition

$$\begin{array}{ccc}\multicolumn{1}{c}{c(t=0,x)=c_0(x)\mathrm{for}x\in \mathrm{\Omega },}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(4)

with $c_0(x)\ge 0$ almost everywhere in Ω.

A fundamental structural property of the systems (1)–(2) is that it represents the constrained $L^2$ -gradient flow associated with the energy functional

$$\mathcal{E}[c]:={\int }_{\mathrm{\Omega }}(c+r)|\mathrm{\nabla }u[c]{|}^2+\frac{\nu }{\gamma }{c}^{\gamma }\mathrm{d}\mathrm{x},$$

(5)

where $u[c]$ is a solution of the Poisson equation (1) subject to homogeneous Neumann boundary condition (3). We note that this property only holds when the source term $S=S(x)$ does not depend on time. The first term in (5) describes the kinetic (pumping) energy necessary for transporting the medium, while the second term is the metabolic expenditure of the network. It is well known that for $\gamma \ge 1$ , the energy functional is convex, see [1,2]. Based on this fact, we construct a unique solution of (1)–(2) as the corresponding gradient flow. We refer to [3,4] for further details.

A variant of the systems (1)–(2) was proposed and studied in Hu and Cai [5,6], where the local permeability of the porous medium is described in terms of a vector field. Detailed mathematical analysis was carried out in a series of papers [3,7–9] and in [10–12], while its various other aspects were studied in [2,13–16]. Another variant with a diagonal permeability tensor was derived formally in Haskovec et al. [17] from the discrete network formation model [16]. The corresponding continuum energy functional was obtained as a limit of a sequence of discrete problems under refinement of the network, with geometry restricted to rectangular networks in two dimensions. In Haskovec et al. [14], the process was carried out rigorously, in terms of the Γ-limit of a sequence of discrete energy functionals with respect to the strong $L^2$ -topology. In Haskovec et al. [1], this was generalized to general symmetric, positive semi-definite conductivity tensors. The continuum energy functional was obtained as a refinement limit of discrete problems on a sequence of triangular meshes. However, the derivation was only formal. Finally, in Alves and Haskovec [18], a graphon limit has been derived rigorously, where (1) is replaced by an integral equation.

Let us note that as long as we are interested in calculating the minimizers of the energy functional (5), i.e., the steady states of the systems (1)–(2), for metabolic exponents $\gamma >1$ , it is sufficient to work with the scalar conductivity c. Indeed, setting $r=0$ and $\nu =1$ for simplicity, the steady-state condition for (2) reads $c^{\gamma -1}=|\mathrm{\nabla }u{|}^2$ . Inserting into (1) gives the p-Laplacian equation

$$\begin{array}{ccc}\multicolumn{1}{c}{-\mathrm{\nabla }\cdot (|\mathrm{\nabla }u{|}^{p-2}\mathrm{\nabla }u)=S,p=\frac{2\gamma }{\gamma -1}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(6)

Note that the exponent $p=\frac{2\gamma }{\gamma -1}>2$ for $\gamma >1$ and $p\to \infty$ when $\gamma \to 1^{+}$ . Now, the steady-state condition for the model [1] with tensor-valued conductivity $C=C(x)\in {\mathrm{ℝ}}^{d\times d}$ reads $|C{|}^{\gamma -2}C=\mathrm{\nabla }u\otimes \mathrm{\nabla }u$ . Therefore, C is the rank-one tensor $C=|\mathrm{\nabla }u{|}^{-2\frac{\gamma -2}{\gamma -1}}\mathrm{\nabla }u\otimes \mathrm{\nabla }u$ , and the flux reads

$$\begin{array}{ccc}\multicolumn{1}{c}{q=-C\mathrm{\nabla }u=-|\mathrm{\nabla }u{|}^{\frac{2}{\gamma -1}}\mathrm{\nabla }u.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Inserting into the Poisson equation $\mathrm{\nabla }\cdot q=S$ then gives (6). Consequently, the equilibrium fluxes of both the tensor- and scalar-valued models coincide and are given by $q=-|\mathrm{\nabla }u{|}^{\frac{2}{\gamma -1}}\mathrm{\nabla }u$ , where u solves (6). We conclude that it is sufficient to work with the scalar-valued conductivity.

Among the most successful numerical solvers for the p-Laplacian are preconditioned gradient descent algorithms [19], barrier methods [20], and quasi-Newton methods [21]. Within this landscape, relaxed Kačanov iterations [22,23] can be considered the gold standard, as they combine robustness with provable convergence. Semi-implicit formulations arising from gradient flows associated with the p-Laplacian energy have also been investigated [24]. In contrast, the gradient flow proposed here serves as a dynamical relaxation: time integration regularizes the nonlinearity and drives the solution to the p-Laplacian steady state, offering an alternative to classical stationary iterations and complementing existing relaxation techniques.

Numerical experiments validate our approach as a competitive alternative to existing relaxation techniques, confirming optimal convergence rates for manufactured solutions and demonstrating robustness for large values of p without the need for adaptive mesh refinement. Numerical evidence indicates that the algorithm can also be applied in the presence of mixed, non-homogeneous Dirichlet and Neumann boundary conditions for (1), and that it exhibits mesh-independent convergence, in the sense that the number of required time steps, nonlinear iterations, and inner linear solver steps remains essentially unaffected by mesh resolution.

This paper is structured as follows. In Section 2, we rigorously derive the energy functional (5) with $\gamma >1$ from the discrete model [5,6] in the refinement Γ-limit of a sequence of triangulations. We also construct its global minimizers in the class of uniformly Lipschitz continuous functions as limits of minimizers of the discrete problem. In Section 3, we consider the numerical approximation of the systems (1)–(2) with general boundary conditions. First, in Section 3.1, we consider the case $\gamma \in (0,1)$ and demonstrate its ability to generate two-dimensional network patterns. Then, in Section 3.2, we solve (1)–(2) with $\gamma >1$ until equilibrium to obtain solutions of the p-Laplacian equation (6).

2. Derivation from the discrete model in two dimensions

The goal of this section is to provide a derivation of the continuum energy functional (5) with $\gamma >1$ as a limit of discrete energy functionals [5,6] posed on a sequence of equilateral triangulations of a bounded polygonal domain $\mathrm{\Omega }\subset {\mathrm{ℝ}}^2$ . Without loss of generality, we will assume that $r>0$ is a constant. For the sake of simplicity, we will consider homogeneous Neumann boundary conditions in the derivation. Nevertheless, the energy functional (5) induces the same gradient-flow structure also for non-homogeneous Neumann boundary conditions and for homogeneous Dirichlet boundary conditions.

We construct a sequence of undirected graphs ${\mathbb{G}}^n=({\mathbb{V}}^n,{\mathbb{E}}^n)$ , $n\in \mathrm{\mathbb{N}}$ , consisting of finite sets of vertices ${\mathbb{V}}^n$ and edges ${\mathbb{E}}^n$ , which are all equilateral triangulations of Ω with edge length $h>0$ approaching zero as $n\to \infty$ . For any $n\in \mathrm{\mathbb{N}}$ , we shall denote ${\mathcal{T}}^n$ the set of all the triangles in the corresponding equilateral triangulation. One possible approach is to iteratively refine each triangle into four equilateral triangles with half the original edge lengths. This would imply $h=\mathcal{O}(2^{-n})$ and the number of triangles (and edges, nodes) proportional to $2^{2n}$ . However, the methods presented in this paper work for any sequence of equilateral triangulations such that $h\to 0$ as $n\to \infty$ . In fact, our approach can be generalized to quasi-uniform triangulations. However, this would significantly increase the technical complexity of the exposition. We therefore choose to restrict to the equilateral case.

Each edge $(i,j)\in {\mathbb{E}}^n$ , connecting the vertices i and $j\in {\mathbb{V}}^n$ , has an associated non-negative conductivity denoted by ${\mathcal{C}}_{\mathrm{ij}}={\mathcal{C}}_{\mathrm{ji}}\ge 0$ . Moreover, each vertex $i\in {\mathbb{V}}^n$ has the pressure $U_i\in \mathrm{ℝ}$ of the material flowing through it. The local mass conservation in each vertex is expressed in terms of the Kirchhoff law

$$\begin{array}{ccc}\multicolumn{1}{c}{-\sum _{j\in \mathcal{N}(i)}({\mathcal{C}}_{\mathrm{ij}}+r)\frac{U_j-U_i}{h}=S_i\mathrm{for}\mathrm{all}i\in {\mathbb{V}}^n,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(7)

where $\mathcal{N}(i)$ denotes the set of vertices connected to $i\in {\mathbb{V}}^n$ through an edge. Moreover, $S={(S_i)}_{i\in {\mathbb{V}}^n}$ is the prescribed vector of strengths of flow sources ( $S_i>0$ ) or sinks ( $S_i<0$ ). Given the matrix of edge conductivities $\mathcal{C}={({\mathcal{C}}_{\mathrm{ij}})}_{(i,j)\in {\mathbb{E}}^n}$ , the Kirchhoff law (7) is a linear system of equations for the vector of pressures $U={(U_i)}_{i\in {\mathbb{V}}^n}$ . A necessary condition for its solvability is the global mass balance

$$\begin{array}{ccc}\multicolumn{1}{c}{\sum _{i\in {\mathbb{V}}^n}{S}_i=0,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(8)

whose validity we assume in the sequel. Note that the solution $U={(U_i)}_{i\in {\mathbb{V}}^n}$ is unique up to an additive constant.

The conductivities ${\mathcal{C}}_{\mathrm{ij}}$ are subject to a minimization process with the energy functional

$$\begin{array}{ccc}\multicolumn{1}{c}{E^n[\mathcal{C}]:=h\sum _{(i,j)\in {\mathbb{E}}^n}({\mathcal{C}}_{\mathrm{ij}}+r)\frac{{(U_j-U_i)}^2}{h^2}+\frac{\nu }{\gamma }{\mathcal{C}}_{\mathrm{ij}}^{\gamma }.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(9)

The first term in the summation describes the (physical) pumping power necessary for transporting the material through the network, while the second term describes the metabolic cost of maintaining the network structure.

For convenience, we shall use both notations $(i,j)\in {\mathbb{E}}^n$ to address edges by the indices of their adjacent vertices, and $e_{\mathrm{ij}}\in {\mathbb{E}}^n$ to address the linear segment in ${\mathrm{ℝ}}^2$ connecting them (i.e., the set of all convex combinations of the vertex coordinates $x_i$ , $x_j\in {\mathrm{ℝ}}^2$ ). For each edge $e_{\mathrm{ij}}\in {\mathbb{E}}^n$ we denote by ${◇}_{\mathrm{ij}}^n\subset {\mathrm{ℝ}}^2$ the adjacent closed diamond, i.e.,

$$\begin{array}{ccc}\multicolumn{1}{c}{{◇}_{\mathrm{ij}}^n:=\left\{x\in \mathrm{\Omega };\mathrm{dist}(x,e_{\mathrm{ij}})\le \mathrm{dist}(x,ẽ)\mathrm{for}\mathrm{all}ẽ\in {\mathbb{E}}^n\right\}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Since the graph ${\mathbb{G}}^n=({\mathbb{V}}^n,{\mathbb{E}}^n)$ represents an equilateral triangulation of Ω, the diamonds are parallelograms with angles $60^{\circ }$ and $120^{\circ }$ . For edges $e_{\mathrm{ij}}$ belonging to the boundary $\partial \mathrm{\Omega }$ the definition of ${◇}_{\mathrm{ij}}$ has to be adjusted accordingly, such that the corresponding diamond is also a parallelogram with angles $60^{\circ }$ and $120^{\circ }$ . For each vertex $i\in {\mathbb{V}}^n$ , we denote ${\mathrm{♢}}_i^n$ the union of adjacent diamonds,

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathrm{♢}}_i^n:=\bigcup _{j\in \mathcal{N}(i)}{◇}_{\mathrm{ij}}^n.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

The main idea of our construction is to map the vector of discrete edge conductivities onto the set of piecewise constant scalar fields on Ω through the mapping ${\mathbb{Q}}^n$ defined as

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathbb{Q}}^n[\mathcal{C}](x):={\mathcal{C}}_{\mathrm{ij}}\mathrm{for}\mathrm{all}x\in {◇}_{\mathrm{ij}}^n.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(10)

${\mathbb{Q}}^n[\mathcal{C}]$ is a constant scalar on each diamond ${◇}_{\mathrm{ij}}^n$ , taking the value ${\mathcal{C}}_{\mathrm{ij}}$ . We also introduce the operator ${\mathbb{Z}}^n$ that maps the edge conductivities onto piecewise constant scalar fields on the triangles $T\in {\mathcal{T}}^n$ ,

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathbb{Z}}^n[\mathcal{C}](x):=\frac{1}{3}\sum _{(i,j)\in T}{\mathcal{C}}_{\mathrm{ij}}\mathrm{for}\mathrm{all}x\in T,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(11)

where the sum goes through the three edges constituting the triangle T. We also introduce the piecewise constant vector field $X^n\in {\mathrm{ℝ}}^2$ , defined as

$$\begin{array}{ccc}\multicolumn{1}{c}{X^n(x):=\frac{x_i-x_j}{|x_i-x_j|}\mathrm{for}\mathrm{all}x\in {◇}_{\mathrm{ij}}^n,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(12)

where we recall that $x_i\in \mathrm{\Omega }$ (respectively $x_j\in \mathrm{\Omega }$ ) denotes the spatial coordinates of the vertices $i\in {\mathbb{V}}^n$ (respectively $j\in {\mathbb{V}}^n$ ).

Our strategy is, by means of the mappings ${\mathbb{Q}}^n$ and ${\mathbb{Z}}^n$ , to establish a connection between a properly rescaled version of the discrete energy functional (9) and its continuum counterpart (5). Following the strategy of [14], the connection will be established through an intermediate functional ${\mathcal{E}}^n$ , where the pressure is a solution of a finite element discretization of the Poisson equation on the triangulation ${\mathcal{T}}^n$ . We shall call this functional the semi-discrete energy functional.

2.1. The Kirchhoff law and the Poisson equation

As the first step, we study the connection between the Kirchhoff law (7) and the Poisson equation (1). We consider the weak formulation of the Poisson equation (1) with the non-negative permeability $c\in L^{\infty }(\mathrm{\Omega })$ and homogeneous Neumann boundary condition,

$$\begin{array}{ccc}\multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}(c+r)\mathrm{\nabla }u\cdot \mathrm{\nabla }\psi \mathrm{d}\mathrm{x}={\int }_{\mathrm{\Omega }}\mathrm{S\psi }\mathrm{d}\mathrm{x}\mathrm{for}\mathrm{all}\psi \in L^2(\mathrm{\Omega }).}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(13)

Since throughout Section 2 we are considering homogeneous Neumann boundary conditions, we impose the global mass balance

$$\begin{array}{ccc}\multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}S(x)\mathrm{d}\mathrm{x}=0.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(14)

We first establish the solvability of (13) with a square-integrable permeability kernel $c\in L^2(\mathrm{\Omega })$ . In the sequel, we denote $H_{\star }^1(\mathrm{\Omega })$ the set of functions in the Sobolev space $H^1(\mathrm{\Omega })$ with zero mean, i.e.

$$\begin{array}{ccc}\multicolumn{1}{c}{H_{\star }^1:=\{u\in H^1(\mathrm{\Omega })\mathrm{s.t.}{\int }_{\mathrm{\Omega }}u\mathrm{d}\mathrm{x}=0\}}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Proposition 1. Let $r>0$ , $S\in L^2(\mathrm{\Omega })$ satisfying the global mass balance (14), and $c\in L^2(\mathrm{\Omega })$ with $c\ge 0$ almost everywhere on Ω. Then, there exists a unique $u\in H_{\star }^1(\mathrm{\Omega })$ verifying (13) for all test functions $\psi \in L^{\infty }(\mathrm{\Omega })$ . Moreover, we have

$$\begin{array}{ccc}\multicolumn{1}{c}{{\parallel \mathrm{\nabla }u\parallel }_{L^2(\mathrm{\Omega })}\le \frac{C_{\mathrm{\Omega }}}{r}{\parallel S\parallel }_{L^2(\mathrm{\Omega })},}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(15)

and

$$\begin{array}{ccc}\multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}(c+r)|\mathrm{\nabla }u{|}^2\mathrm{d}\mathrm{x}={\int }_{\mathrm{\Omega }}\mathrm{Su}\mathrm{d}\mathrm{x}\le \frac{C_{\mathrm{\Omega }}^2}{r}{\parallel S\parallel }_{L^2(\mathrm{\Omega })}^2,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(16)

where $C_{\mathrm{\Omega }}$ is the Poincaré constant of the domain Ω.

The proof is obtained as a slight modification of the proof of [3, Lemma 6] or [18, Lemma 5.5]. Now, for the non-negative permeability $c\in L^{\infty }(\mathrm{\Omega })$ , we consider the modified Poisson equation

$$\begin{array}{ccc}\multicolumn{1}{c}{-2\mathrm{\nabla }\cdot ((c+r)(X^n\otimes X^n)\mathrm{\nabla }u)=S,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(17)

with $X^n$ defined in (12) and the symbol⊗ denoting the standard tensor product. The factor 2 is introduced because, as we shall see, $X^n\otimes X^n\rightharpoonup \frac{1}{2}\mathbb{I}$ as $n\to \infty$ . In the sequel, we shall work with the first-order $H^1$ finite element approximation of (17),

$$\begin{array}{ccc}\multicolumn{1}{c}{2{\int }_{\mathrm{\Omega }}(c+r)\mathrm{\nabla }\psi \cdot (X^n\otimes X^n)\mathrm{\nabla }u\mathrm{d}\mathrm{x}={\int }_{\mathrm{\Omega }}\mathrm{S\psi }\mathrm{d}\mathrm{x}\mathrm{for}\mathrm{all}\psi \in Y_1^n.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(18)

here $Y_1^n$ denotes the space of continuous, piecewise linear functions on ${\mathcal{T}}^n$ ,

$$\begin{array}{ccc}\multicolumn{1}{c}{Y_1^n:=\left\{u\in C^0(\mathrm{\Omega });u\mathrm{linear}\mathrm{on}\mathrm{each}T\in {\mathcal{T}}^n\right\}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Moreover, we denote $Y_0^n$ the space of bounded, piecewise constant functions on the triangulation ${\mathcal{T}}^n$ ,

$$\begin{array}{ccc}\multicolumn{1}{c}{Y_0^n:=\left\{u\in L^{\infty }(\mathrm{\Omega });u\mathrm{constant}\mathrm{on}\mathrm{each}T\in {\mathcal{T}}^n\right\}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Existence of solutions of (18) is established by a standard application of the Lax–Milgram theorem, where coercivity of the respective bilinear form with $r>0$ follows from formula (22) below. The following proposition establishes a connection between the solution of (18) with $c:={\mathbb{Q}}^n[\mathcal{C}]$ and the Kirchhoff law (7).

Proposition 2. Let $S\in L^2(\mathrm{\Omega })$ verify the global mass balance assumption (14) and fix $n\in \mathrm{\mathbb{N}}$ . Let ${\mathcal{C}}_{\mathrm{ij}}\ge 0$ for all $(i,j)\in {\mathbb{E}}^n$ and denote $c(x):={\mathbb{Q}}^n[\mathcal{C}](x)$ for $x\in \mathrm{\Omega }$ . Let $u^n:=u[c]\in Y_1^n$ be the solution of the discrete Poisson equation (18) with conductivity $c=c(x)$ . For $i\in {\mathbb{V}}^n$ denote $U_i^n:=u^n(x_i)$ .

Then, we have, for all $i\in {\mathbb{V}}^n$ ,

$$\begin{array}{ccc}\multicolumn{1}{c}{-\sum _{j\in \mathcal{N}(i)}({\mathcal{C}}_{\mathrm{ij}}+r)\frac{U_j^n-U_i^n}{h}=S_i^n,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(19)

where we denoted

$$\begin{array}{ccc}\multicolumn{1}{c}{S_i^n:=\frac{h}{2\mathrm{vol}({◇}_h)}{\int }_{\mathrm{\Omega }}S{\psi }_i^n\mathrm{d}\mathrm{x}\mathrm{for}\mathrm{all}i\in {\mathbb{V}}^n,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(20)

and ${\psi }_i^n\in Y_1^n$ is a piecewise linear function supported on ${\mathrm{♢}}_i$ , taking the vertex values ${\psi }_i^n(x_i)=1$ and ${\psi }_i^n(x_j)=0$ for all $j\ne i$ .

The proof is a minor modification of the proof of [1, Proposition 2.1]. Note that (19) is the Kirchhoff law (7) with right-hand side given by (20). We have

$$\begin{array}{ccccccc}\multicolumn{1}{c}{\frac{2\mathrm{vol}({◇}_h)}{h}\sum _{i\in {\mathbb{V}}^n}{S}_i^n}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{\sum _{i\in {\mathbb{V}}^n}{\int }_{\mathrm{\Omega }}S(x){\psi }_i^n(x)\mathrm{d}\mathrm{x}}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{\sum _{T\in {\mathcal{T}}^n}{\int }_{T}S(x)\sum _{i\in {\mathbb{V}}^n}{\psi }_i^n(x)\mathrm{d}\mathrm{x}}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}S(x)\mathrm{d}\mathrm{x},}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

where we used the fact that the sum of barycentric coordinates on each triangle is constant unity, i.e., ${\mathrm{\sum }}_{i\in {\mathbb{V}}^n}{\psi }_i^n(x)\equiv 1$ for all $x\in \mathrm{\Omega }$ . Consequently, assumption (14) implies the discrete mass balance (8).

Moreover, for equilateral triangulations, we have

$$\begin{array}{ccc}\multicolumn{1}{c}{\frac{h}{2\mathrm{vol}({◇}_h)}=\frac{\sqrt{3}}{h},{\int }_{\mathrm{\Omega }}{({\psi }_i^n)}^2\mathrm{d}\mathrm{x}=\frac{\sqrt{3}}{24}{h}^2.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Consequently, by the Cauchy–Schwarz inequality,

$$\begin{array}{ccc}\multicolumn{1}{c}{{(S_i^n)}^2\le \frac{\sqrt{3}}{8}{\int }_{{\mathrm{♢}}_i}S{(x)}^2\mathrm{d}\mathrm{x},}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

so that

$$\begin{array}{ccc}\multicolumn{1}{c}{\sum _{i\in {\mathbb{V}}^n}{(S_i^n)}^2\le \frac{3\sqrt{3}}{8}{\int }_{\mathrm{\Omega }}S{(x)}^2\mathrm{d}\mathrm{x}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(21)

To facilitate the limit passage $n\to \infty$ , we shall work with the discrete Poisson equation (18) with $c:={\mathbb{Z}}^n[\mathcal{C}]$ . We have the following result.

Lemma 3. Fix $n\in \mathrm{\mathbb{N}}$ . For any $u,v\in Y_1^n(\mathrm{\Omega })$ and any equilateral triangle $T\in {\mathcal{T}}^n$ , we have

$$\begin{array}{ccc}\multicolumn{1}{c}{{\int }_T\mathrm{\nabla }u\cdot (X^n\otimes X^n)\mathrm{\nabla }v\mathrm{d}\mathrm{x}=\frac{1}{2}{\int }_T\mathrm{\nabla }u\cdot \mathrm{\nabla }v\mathrm{d}\mathrm{x}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(22)

Proof: Let us denote the three edges of the triangle T by $e_1$ , $e_2$ , $e_3$ and, without loss of generality, assume $|e_1|=|e_2|=|e_3|=1$ . Since both $\mathrm{\nabla }u$ and $\mathrm{\nabla }v$ are, by definition, constant on the triangle, with (12), we have

$$\begin{array}{ccccccc}\multicolumn{1}{c}{{\int }_T\mathrm{\nabla }u\cdot (X^n\otimes X^n)\mathrm{\nabla }v\mathrm{d}\mathrm{x}}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{\mathrm{\nabla }u\cdot \left({\int }_T{X}^n\otimes X^n\mathrm{d}\mathrm{x}\right)\mathrm{\nabla }v}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{\mathrm{\nabla }u\cdot \left(\frac{1}{4\sqrt{3}}\sum _{i=1}^3{e}_i\otimes e_i\right)\mathrm{\nabla }v.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

By rotational symmetry, the matrix ${\mathrm{\sum }}_{i=1}^3{e}_i\otimes e_i$ is a scalar multiple of the identity. Since

$$\begin{array}{ccc}\multicolumn{1}{c}{\mathrm{trace}\left(\sum _{i=1}^3{e}_i\otimes e_i\right)=\sum _{i=1}^3|e_i{|}^2=3,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

we conclude ${\mathrm{\sum }}_{i=1}^3{e}_i\otimes e_i=\frac{3}{2}\mathbb{I}$ . Consequently,

$$\begin{array}{ccc}\multicolumn{1}{c}{{\int }_T\mathrm{\nabla }u\cdot (X^n\otimes X^n)\mathrm{\nabla }v\mathrm{d}\mathrm{x}=\frac{\sqrt{3}}{8}\mathrm{\nabla }u\cdot \mathrm{\nabla }v=\frac{1}{2}{\int }_T\mathrm{\nabla }u\cdot \mathrm{\nabla }v\mathrm{d}\mathrm{x}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

■

Lemma 3 allows us to write the left-hand side of (18) with $c:={\mathbb{Z}}^n[\mathcal{C}]$ as

$$\begin{array}{ccccccc}\multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}({\mathbb{Z}}^n[\mathcal{C}]+r)\mathrm{\nabla }\psi \cdot (X^n\otimes X^n)\mathrm{\nabla }u\mathrm{d}\mathrm{x}}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{\sum _{T\in {\mathcal{T}}^n}({\mathbb{Z}}^n{[\mathcal{C}]}_T+r){\int }_T\mathrm{\nabla }\psi \cdot (X^n\otimes X^n)\mathrm{\nabla }u\mathrm{d}\mathrm{x}}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{\frac{1}{2}{\int }_{\mathrm{\Omega }}({\mathbb{Z}}^n[\mathcal{C}]+r)\mathrm{\nabla }\psi \cdot \mathrm{\nabla }u\mathrm{d}\mathrm{x},}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

where ${\mathbb{Z}}^n{[\mathcal{C}]}_T$ denotes the constant value of ${\mathbb{Z}}^n[\mathcal{C}]$ on the triangle $T\in {\mathcal{T}}^n$ . Consequently, (18) with $c:={\mathbb{Z}}^n[\mathcal{C}]$ can be equivalently formulated as

$$\begin{array}{ccc}\multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}({\mathbb{Z}}^n[\mathcal{C}]+r)\mathrm{\nabla }\psi \cdot \mathrm{\nabla }u\mathrm{d}\mathrm{x}={\int }_{\mathrm{\Omega }}\mathrm{S\psi }\mathrm{d}\mathrm{x}\mathrm{for}\mathrm{all}\psi \in Y_1^n.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(23)

This formulation facilitates the passage to the Γ-limit as $n\to +\infty$ since it avoids the oscillatory term $X^n\otimes X^n$ that converges only weakly to $\frac{1}{2}\mathbb{I}$ .

To estimate the difference between the scalar fields ${\mathbb{Q}}^n[\mathcal{C}]$ and ${\mathbb{Z}}^n[\mathcal{C}]$ , let us for any $T\in {\mathcal{T}}^n$ denote

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathbb{D}}^n[\mathcal{C}](T):=\max \{|{\mathcal{C}}_{\mathrm{ij}}-{\mathcal{C}}_{\mathrm{ik}}|,|{\mathcal{C}}_{\mathrm{ik}}-{\mathcal{C}}_{\mathrm{jk}}|,|{\mathcal{C}}_{\mathrm{ij}}-{\mathcal{C}}_{\mathrm{jk}}|\}}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(24)

where $(i,j,k)$ are the vertices constituting T.

Lemma 4. For any $\mathcal{C}\in {\mathrm{ℝ}}^{|{\mathbb{E}}^n|}$ , we have

$$\begin{array}{ccc}\multicolumn{1}{c}{{\parallel {\mathbb{Q}}^n[\mathcal{C}]-{\mathbb{Z}}^n[\mathcal{C}]\parallel }_{L^{\infty }(\mathrm{\Omega })}\le \frac{2}{3}\underset{T\in {\mathcal{T}}^n}{\max }{\mathbb{D}}^n[\mathcal{C}](T).}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(25)

Proof: Let us fix any triangle $T\in {\mathcal{T}}^n$ and denote its vertices i, j, $k\in {\mathbb{V}}^n$ . For the half-diamond ${◇}_{\mathrm{ij}}\cap T$ , we have

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathbb{Q}}^n[\mathcal{C}]\equiv {\mathcal{C}}_{\mathrm{ij}},{\mathbb{Z}}^n[\mathcal{C}]\equiv \frac{1}{3}({\mathcal{C}}_{\mathrm{ij}}+{\mathcal{C}}_{\mathrm{jk}}+{\mathcal{C}}_{\mathrm{ik}}).}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Consequently,

$$\begin{array}{ccc}\multicolumn{1}{c}{|{\mathbb{Q}}^n[{\mathcal{C}}^n]-{\mathbb{Z}}^n[{\mathcal{C}}^n]|=\frac{1}{3}|2{\mathcal{C}}_{\mathrm{ij}}-{\mathcal{C}}_{\mathrm{jk}}-{\mathcal{C}}_{\mathrm{ik}}|\le \frac{2}{3}\max \left\{|{\mathcal{C}}_{\mathrm{ij}}-{\mathcal{C}}_{\mathrm{jk}}|,|{\mathcal{C}}_{\mathrm{ij}}-{\mathcal{C}}_{\mathrm{ik}}|\right\}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Repeating the estimate for the other two half-diamonds constituting T, and taking the maximum over all $T\in {\mathcal{T}}^n$ , we obtain (25). ■

2.2. The semi-discrete energy functionals

For the non-negative permeability $c\in L^{\infty }(\mathrm{\Omega })$ , we introduce the semi-discrete energy functional

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathcal{E}}^n[c]={\int }_{\mathrm{\Omega }}2\left(c+r\right){|X^n\cdot \mathrm{\nabla }u[c]|}^2+\frac{\nu }{\gamma }{c}^{\gamma }\mathrm{d}\mathrm{x},}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(26)

where $u[c]\in Y_1^n$ the unique weak solution of the discrete Poisson equation (18). Let us note that using the first-order $H^1$ finite element approximation (18) of the modified Poisson equation (17) in the definition of the semi-discrete energy functional (26) is a fundamental step in the derivation of the continuum limit because it allows us to establish a connection to the discrete model. Namely, we introduce the rescaled version of the discrete energy functional (9),

$$\begin{array}{ccc}\multicolumn{1}{c}{{\overline{E}}^n[\mathcal{C}]=\mathrm{vol}({◇}_h)\sum _{(i,j)\in {\mathbb{E}}^n}2\left({\mathcal{C}}_{\mathrm{ij}}+r\right)\frac{{(U_i-U_j)}^2}{h^2}+\frac{\nu }{\gamma }{\mathcal{C}}_{\mathrm{ij}}^{\gamma },}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(27)

where $U=U[\mathcal{C}]$ is a solution of the rescaled Kirchhoff law (19). The following proposition, which is a slight modification of [1, Proposition 2.2], establishes the equality of ${\mathcal{E}}^n$ evaluated at the piecewise constant function ${\mathbb{Q}}^n[\mathcal{C}]$ to the value of ${\overline{E}}^n[\mathcal{C}]$ .

Proposition 5. For any $n\in \mathrm{\mathbb{N}}$ and any vector of non-negative conductivities $\mathcal{C}={({\mathcal{C}}_{\mathrm{ij}})}_{(i,j)\in {\mathbb{E}}^n}$ , we have

$$\begin{array}{ccc}\multicolumn{1}{c}{{\overline{E}}^n[\mathcal{C}]={\mathcal{E}}^n\left[{\mathbb{Q}}^n[\mathcal{C}]\right],}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(28)

with the rescaled discrete energy functional ${\overline{E}}^n$ defined in (27) and the semi-discrete energy ${\mathcal{E}}^n$ given by (26).

We now estimate the difference of energies calculated with $c:={\mathbb{Q}}^n[\mathcal{C}]$ and $c:={\mathbb{Z}}^n[\mathcal{C}]$ .

Lemma 6. Fix $n\in \mathrm{\mathbb{N}}$ . For any $\mathcal{C}\in {\mathrm{ℝ}}^{|{\mathbb{E}}^n|}$ , we have

$$\begin{array}{c}\multicolumn{1}{c}{|{\mathcal{E}}^n[{\mathbb{Q}}^n[\mathcal{C}]]-{\mathcal{E}}^n[{\mathbb{Z}}^n[\mathcal{C}]]|}\\ \multicolumn{1}{c}{\le \frac{2}{3}\underset{T\in {\mathcal{T}}^n}{\max }{\mathbb{D}}^n[\mathcal{C}](T)\left(2{\parallel \mathrm{\nabla }{u}^{\mathbb{Q}}\parallel }_{L^2(\mathrm{\Omega })}{\parallel \mathrm{\nabla }{u}^{\mathbb{Z}}\parallel }_{L^2(\mathrm{\Omega })}+|\mathrm{\Omega }|\nu {\parallel \mathcal{C}\parallel }_{{\ell }_{\infty }}^{\gamma -1}\right),}\end{array}$$

(29)

where we denoted $u^{\mathbb{Q}}\in Y_1^n$ the solution of the discrete Poisson equation (18) with permeability ${\mathbb{Q}}^n[\mathcal{C}]$ , and $u^{\mathbb{Z}}\in Y_1^n$ the solution of (18) with permeability ${\mathbb{Z}}^n[\mathcal{C}]$ .

Proof: Using $u^{\mathbb{Q}}$ as a test function in the discrete Poisson equation (18) for $u^{\mathbb{Q}}$ yields for the kinetic part of (26),

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathcal{E}}_{\mathrm{kin}}^{\mathrm{n}}[{\mathbb{Q}}^n[\mathcal{C}]]=2{\int }_{\mathrm{\Omega }}{\mathbb{Q}}^n[\mathcal{C}]|X^n\cdot \mathrm{\nabla }{u}^{\mathbb{Q}}{|}^2\mathrm{d}\mathrm{x}=\int S{u}^{\mathbb{Q}}\mathrm{d}\mathrm{x}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Analogously, we have

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathcal{E}}_{\mathrm{kin}}^{\mathrm{n}}[{\mathbb{Z}}^n[\mathcal{C}]]=2{\int }_{\mathrm{\Omega }}{\mathbb{Z}}^n[\mathcal{C}]|X^n\cdot \mathrm{\nabla }{u}^{\mathbb{Z}}{|}^2\mathrm{d}\mathrm{x}=\int S{u}^{\mathbb{Z}}\mathrm{d}\mathrm{x}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Therefore,

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathcal{E}}_{\mathrm{kin}}^{\mathrm{n}}[{\mathbb{Q}}^n[\mathcal{C}]]-{\mathcal{E}}_{\mathrm{kin}}^{\mathrm{n}}[{\mathbb{Z}}^n[\mathcal{C}]]={\int }_{\mathrm{\Omega }}S\left(u^{\mathbb{Q}}-u^{\mathbb{Z}}\right)\mathrm{d}\mathrm{x}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Now, using $u^{\mathbb{Z}}$ as a test function in the discrete Poisson equation for $u^{\mathbb{Q}}$ and vice versa, and subtracting the resulting identities, we obtain

$$\begin{array}{ccc}\multicolumn{1}{c}{2{\int }_{\mathrm{\Omega }}\left({\mathbb{Q}}^n[\mathcal{C}]-{\mathbb{Z}}^n[\mathcal{C}]\right)(X^n\cdot \mathrm{\nabla }{u}^{\mathbb{Q}})(X^n\cdot \mathrm{\nabla }{u}^{\mathbb{Z}})\mathrm{d}\mathrm{x}={\int }_{\mathrm{\Omega }}S\left(u^{\mathbb{Q}}-u^{\mathbb{Z}}\right)\mathrm{d}\mathrm{x}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

We combine the last two identities and use the Cauchy–Schwarz inequality to estimate

$$\begin{array}{ccccccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{\left|{\int }_{\mathrm{\Omega }}\left({\mathbb{Q}}^n[\mathcal{C}]-{\mathbb{Z}}^n[\mathcal{C}]\right)(X^n\cdot \mathrm{\nabla }{u}^{\mathbb{Q}})(X^n\cdot \mathrm{\nabla }{u}^{\mathbb{Z}})\mathrm{d}\mathrm{x}\right|}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{\le {\parallel {\mathbb{Q}}^n[\mathcal{C}]-{\mathbb{Z}}^n[\mathcal{C}]\parallel }_{L^{\infty }(\mathrm{\Omega })}{\parallel \mathrm{\nabla }{u}^{\mathbb{Q}}\parallel }_{L^2(\mathrm{\Omega })}{\parallel \mathrm{\nabla }{u}^{\mathbb{Z}}\parallel }_{L^2(\mathrm{\Omega })},}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

where we also used the uniform bound $|X^n|\le 1$ .

To estimate the difference of the metabolic parts of (26), we use the inequality $|a^{\gamma }-b^{\gamma }|\le \gamma \max {\{a,b\}}^{\gamma -1}|a-b|$ , which holds for a, $b\ge 0$ and $\gamma \ge 1$ . We obtain

$$\begin{array}{ccc}\multicolumn{1}{c}{\frac{\nu }{\gamma }{\int }_{\mathrm{\Omega }}|{\mathbb{Q}}^n{[\mathcal{C}]}^{\gamma }-{\mathbb{Z}}^n{[\mathcal{C}]}^{\gamma }|\mathrm{d}\mathrm{x}\le |\mathrm{\Omega }|\nu {\parallel \mathcal{C}\parallel }_{{\ell }_{\infty }}^{\gamma -1}{\parallel {\mathbb{Q}}^n[\mathcal{C}]-{\mathbb{Z}}^n[\mathcal{C}]\parallel }_{L^{\infty }(\mathrm{\Omega })}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

We conclude by using (25). ■

With the uniform $H^1$ -bound (15) on the solutions of the Poisson equation (18) provided by Proposition 1, we readily obtain the following corollary.

Corollary 7. Let ${\mathcal{C}}^n\in {\mathrm{ℝ}}^{|{\mathbb{E}}^n|}$ , $n\in \mathrm{\mathbb{N}}$ , be a sequence of non-negative conductivities such that

$$\begin{array}{ccc}\multicolumn{1}{c}{\underset{n\in \mathrm{\mathbb{N}}}{\sup }{\parallel {\mathcal{C}}^n\parallel }_{{\ell }_{\infty }}<\infty \mathrm{and}\underset{n\to \infty }{\lim }\underset{T\in {\mathcal{T}}^n}{\max }{\mathbb{D}}^n[{\mathcal{C}}^n](T)=0.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(30)

Then

$$\begin{array}{ccc}\multicolumn{1}{c}{\underset{n\to \infty }{\lim }|{\mathcal{E}}^n[{\mathbb{Q}}^n[{\mathcal{C}}^n]]-{\mathcal{E}}^n[{\mathbb{Z}}^n[{\mathcal{C}}^n]]|=0.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

2.3. The continuum energy functional

In this section, we study the properties of the continuum energy functional (5) and its relation to its semi-discrete version (26). We start by proving the continuity of (5) with respect to the norm topology of $L^{\omega }(\mathrm{\Omega })$ , $\omega :=\max \{2,\gamma \}$ .

Lemma 8. Let $S\in L^2(\mathrm{\Omega })$ satisfy (14) and fix $r>0$ . Let the sequence ${(c^n)}_{n\in \mathrm{\mathbb{N}}}\subset L^{\omega }(\mathrm{\Omega })$ converge to $c\in L^{\omega }(\mathrm{\Omega })$ in the norm topology of $L^{\omega }(\mathrm{\Omega })$ as $n\to \infty$ . Then

$$\begin{array}{ccc}\multicolumn{1}{c}{\mathcal{E}[c]=\underset{n\to \infty }{\lim }\mathcal{E}[c^n].}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Proof: Let ${(u^n)}_{n\in \mathrm{\mathbb{N}}}\subset L_0^2(\mathrm{\Omega })$ be a sequence of weak solutions of the Poisson equation (13) with permeability $c^n\in L^{\omega }(\mathrm{\Omega })$ , provided by Proposition 1. Identity (16) gives

$$\begin{array}{ccc}\multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}(c^n+r)|\mathrm{\nabla }{u}^n{|}^2\mathrm{d}\mathrm{x}={\int }_{\mathrm{\Omega }}S{u}^n\mathrm{d}\mathrm{x}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(31)

Due to the uniform bound (15), there exists $u\in H_{\star }^1(\mathrm{\Omega })$ , a weak limit of (a subsequence of) $u^n$ , verifying (13) for all test functions $\phi \in L_0^{\infty }(\mathrm{\Omega })$ . Since $c^n\to c$ strongly in $L^2(\mathrm{\Omega })$ , we can pass to the limit in the weak formulation (13) which establishes u as its weak solution with permeability c and test functions $\phi \in L^{\infty }(\mathrm{\Omega })$ . Again, (16) gives

$$\begin{array}{ccc}\multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}(c+r)|\mathrm{\nabla }u{|}^2\mathrm{d}\mathrm{x}={\int }_{\mathrm{\Omega }}\mathrm{Su}\mathrm{d}\mathrm{x}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Combining with (31) we have

$$\begin{array}{ccc}\multicolumn{1}{c}{\underset{n\to \infty }{\lim }{\int }_{\mathrm{\Omega }}(c^n+r)|\mathrm{\nabla }{u}^n{|}^2\mathrm{d}\mathrm{x}=\underset{n\to \infty }{\lim }{\int }_{\mathrm{\Omega }}S{u}^n\mathrm{d}\mathrm{x}={\int }_{\mathrm{\Omega }}\mathrm{Su}\mathrm{d}\mathrm{x}={\int }_{\mathrm{\Omega }}(c+r)|\mathrm{\nabla }u{|}^2\mathrm{d}\mathrm{x},}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

which proves the strong continuity of the kinetic part of the energy functional E. The metabolic part is a multiple of the $L^{\gamma }$ -norm of c; consequently, its strong continuity is obvious. ■

A major observation is that for $\gamma \ge 1$ , the energy functional $\mathcal{E}$ given by (5) is convex. A proof can be found in [1, Proposition 3.2] or [18, Lemma 6.4]. Then, since $\mathcal{E}$ is strongly continuous on $L^{\omega }(\mathrm{\Omega })$ due to Lemma 8, an application of Mazur’s lemma [25] establishes its weak lower semi-continuity, i.e.,

$$\begin{array}{ccc}\multicolumn{1}{c}{\mathcal{E}[c]\le \underset{n\to \infty }{\lim \inf }\mathcal{E}[c^n]}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(32)

whenever ${(c^n)}_{n\in \mathrm{\mathbb{N}}}\subset L^{\omega }(\mathrm{\Omega })$ converges weakly in $L^{\omega }(\mathrm{\Omega })$ to c.

Next, we establish an inequality between the semi-discrete and continuum energy functionals restricted to piecewise constant permeabilities on the triangles $T\in {\mathcal{T}}^n$ .

Lemma 9. For any $n\in \mathrm{\mathbb{N}}$ and non-negative $c^n\in Y_0^n$ , we have

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathcal{E}}^n[c^n]\le \mathcal{E}[c^n],}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(33)

with the semi-discrete energy functional ${\mathcal{E}}^n$ given by (26) and $\mathcal{E}$ given by (5).

Proof: Let us recall the weak formulation of the discrete Poisson equation (18) with test function $\psi \in Y_1^n$ ,

$$\begin{array}{ccc}\multicolumn{1}{c}{2{\int }_{\mathrm{\Omega }}(c^n+r)\mathrm{\nabla }\psi \cdot (X^n\otimes X^n)\mathrm{\nabla }{u}^n\mathrm{d}\mathrm{x}={\int }_{\mathrm{\Omega }}\mathrm{S\psi }\mathrm{d}\mathrm{x}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Lemma 3 gives

$$\begin{array}{ccccccc}\multicolumn{1}{c}{2{\int }_{\mathrm{\Omega }}(c^n+r)\mathrm{\nabla }\psi \cdot (X^n\otimes X^n)\mathrm{\nabla }{u}^n\mathrm{d}\mathrm{x}}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{2\sum _{T\in {\mathcal{T}}^n}\left(c_T^n+r\right){\int }_T(X^n\cdot \mathrm{\nabla }{u}^n)(X^n\cdot \mathrm{\nabla }\psi )\mathrm{d}\mathrm{x}}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{\sum _{T\in {\mathcal{T}}^n}|T|\left(c_T^n+r\right){(\mathrm{\nabla }{u}^n)}_T\cdot {(\mathrm{\nabla }\psi )}_T,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

where $c_T^n$ denotes the constant value of $c^n$ on the triangle $T\in {\mathcal{T}}^n$ , and similarly for ${(\mathrm{\nabla }\psi )}_T$ and ${(\mathrm{\nabla }{u}^n)}_T$ .

Let $u\in H_{\star }^1(\mathrm{\Omega })$ be the unique weak solution of the Poisson equation (1) with permeability $c^n$ . Using again $\psi \in Y_1^n$ as a test function gives

$$\begin{array}{ccc}\multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}(c^n+r)\mathrm{\nabla }u\cdot \mathrm{\nabla }\psi \mathrm{d}\mathrm{x}=\sum _{T\in {\mathcal{T}}^n}\left(c_T^n+r\right){(\mathrm{\nabla }\psi )}_T\cdot {\int }_T\mathrm{\nabla }u\mathrm{d}\mathrm{x}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Due to the uniqueness (in the respective spaces of functions with vanishing mean) of the solutions of (18) and (1), we have

$$\begin{array}{ccc}\multicolumn{1}{c}{{(\mathrm{\nabla }{u}^n)}_T=\frac{1}{|T|}{\int }_T\mathrm{\nabla }u\mathrm{d}\mathrm{x}\mathrm{for}\mathrm{all}T\in {\mathcal{T}}^n.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Then, we have for the kinetic parts of the semi-discrete and continuum energies

$$\begin{array}{ccccccc}\multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}(c^n+r)|\mathrm{\nabla }{u}^n{|}^2\mathrm{d}\mathrm{x}}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{\sum _{T\in {\mathcal{T}}^n}|T|(c_T^n+r)|{(\mathrm{\nabla }{u}^n)}_T{|}^2}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{\sum _{T\in {\mathcal{T}}^n}|T|(c_T^n+r){|\frac{1}{|T|}{\int }_T\mathrm{\nabla }u\mathrm{d}\mathrm{x}|}^2}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{\le }& \multicolumn{1}{c}{\sum _{T\in {\mathcal{T}}^n}(c_T^n+r){\int }_T|\mathrm{\nabla }u{|}^2\mathrm{d}\mathrm{x}}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}(c^n+r)|\mathrm{\nabla }u{|}^2\mathrm{d}\mathrm{x},}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

where we used the Jensen inequality. Since the metabolic part $\frac{\nu }{\gamma }{\int }_{\mathrm{\Omega }}{(c^n)}^{\gamma }\mathrm{d}\mathrm{x}$ is identical for both the semi-discrete and continuum energies, we conclude (33). ■

Lemma 10. Let ${(c^n)}_{n\in \mathrm{\mathbb{N}}}$ , be a sequence of piecewise constant non-negative conductivities $c^n\in Y_0^n$ such that $\mathrm{su}{p}_{n\in \mathrm{\mathbb{N}}}{\parallel c^n\parallel }_{L^{\infty }(\mathrm{\Omega })}<\infty$ . Then

$$\begin{array}{ccc}\multicolumn{1}{c}{\underset{n\to \infty }{\lim }|{\mathcal{E}}^n[c^n]-\mathcal{E}[c^n]|=0.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(34)

Proof: Using the solution $u^n\in H_{\star }^1(\mathrm{\Omega })$ of (13) as a test function, we obtain

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathcal{E}}_{\mathrm{kin}}[c^n]={\int }_{\mathrm{\Omega }}(c^n+r)|\mathrm{\nabla }{u}^n{|}^2\mathrm{d}\mathrm{x}={\int }_{\mathrm{\Omega }}S{u}^n\mathrm{d}\mathrm{x}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

On the contrary, using the solution $u_h^n\in Y_1^n(\mathrm{\Omega })$ of (18) as a test function, we obtain

$$\begin{array}{ccc}\multicolumn{1}{c}{{\mathcal{E}}_{\mathrm{kin}}^{\mathrm{n}}[c^n]={\int }_{\mathrm{\Omega }}2(c^n+r)|X^n\cdot \mathrm{\nabla }{u}_h^n{|}^2\mathrm{d}\mathrm{x}={\int }_{\mathrm{\Omega }}S{u}_h^n\mathrm{d}\mathrm{x}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Consequently,

$$\begin{array}{ccccccc}\multicolumn{1}{c}{\left|{\mathcal{E}}^n[c^n]-\mathcal{E}[c^n]\right|}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{\left|{\int }_{\mathrm{\Omega }}S\left(u^n-u_h^n\right)\mathrm{d}\mathrm{x}\right|}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{\le }& \multicolumn{1}{c}{{\parallel S\parallel }_{L^2(\mathrm{\Omega })}{\parallel u^n-P^n[u^n]\parallel }_{L^2(\mathrm{\Omega })}+|{\int }_{\mathrm{\Omega }}S\left(P^n[u^n]-u_h^n\right)\mathrm{d}\mathrm{x}|,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

where $P^n:H^1(\mathrm{\Omega })\to Y_1^n$ is an $H^1$ -stable interpolator, i.e., Scott–Zhang or Clement. Note that due to the uniform boundedness of $u^n$ in $H_{\star }^1(\mathrm{\Omega })$ provided by (15), we have

$$\begin{array}{ccc}\multicolumn{1}{c}{\underset{n\to \infty }{\lim }{\parallel u^n-P^n[u^n]\parallel }_{H^1(\mathrm{\Omega })}=0.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(35)

Moreover, using $P^n[u^n]\in Y_1^n$ as a test function in the discrete Poisson equation (18) for $u_h^n$ , we have

$$\begin{array}{ccccccc}\multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}2(c^n+r)(X^n\cdot \mathrm{\nabla }{u}_h^n)(X^n\cdot \mathrm{\nabla }{P}^n[u^n])\mathrm{d}\mathrm{x}}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}(c^n+r)\mathrm{\nabla }{u}_h^n\cdot \mathrm{\nabla }{P}^n[u^n]\mathrm{d}\mathrm{x}}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}S{P}^n[u^n]\mathrm{d}\mathrm{x},}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

where we used Lemma 3 in the first line, recalling that $c^n\in Y_0^n$ . On the other hand, using $u_h^n\in Y_1^n$ as a test function in the Poisson equation (13) for $u^n$ , we obtain

$$\begin{array}{ccc}\multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}(c^n+r)\mathrm{\nabla }{u}^n\cdot \mathrm{\nabla }{u}_h^n\mathrm{d}\mathrm{x}={\int }_{\mathrm{\Omega }}S{u}_h^n\mathrm{d}\mathrm{x}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Consequently,

$$\begin{array}{ccccccc}\multicolumn{1}{c}{|{\int }_{\mathrm{\Omega }}S\left(P^n[u^n]-u_h^n\right)\mathrm{d}\mathrm{x}|}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{|{\int }_{\mathrm{\Omega }}(c^n+r)\mathrm{\nabla }{u}_h^n\cdot \mathrm{\nabla }\left(P^n[u^n]-u^n\right)\mathrm{d}\mathrm{x}|}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{\le }& \multicolumn{1}{c}{\underset{n\in \mathrm{\mathbb{N}}}{\sup }{\parallel c^n+r\parallel }_{L^{\infty }(\mathrm{\Omega })}{\parallel \mathrm{\nabla }{u}_h^n\parallel }_{L^2(\mathrm{\Omega })}{\parallel \mathrm{\nabla }\left(P^n[u^n]-u^n\right)\parallel }_{L^2(\mathrm{\Omega })}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Using again the uniform boundedness of $u_h^n$ in $H_{\star }^1(\mathrm{\Omega })$ provided by (15) of Proposition 1, combined with the $H^1$ -stability (35), we conclude. ■

2.4. Γ-convergence of the semi-discrete energy functionals and construction of global energy minimizers

The main result of this section is the Γ-convergence of the sequence of the semi-discrete energy functionals (26) in the appropriate topology.

Theorem 11. Let ${(c^n)}_{n\in \mathrm{\mathbb{N}}}$ be a sequence of non-negative conductivities $c^n\in Y_0^n$ such that $\mathrm{su}{p}_{n\in \mathrm{\mathbb{N}}}{\parallel c^n\parallel }_{L^{\infty }(\mathrm{\Omega })}$ $<\infty$ and $c^n\rightharpoonup c$ weakly- ∗ in $L^{\infty }(\mathrm{\Omega })$ . Then

$$\begin{array}{ccc}\multicolumn{1}{c}{\mathcal{E}[c]\le \underset{n\to \infty }{\lim \inf }{\mathcal{E}}^n[c^n].}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(36)

Moreover, for any $c\in L^{\infty }(\mathrm{\Omega })$ , there exists a sequence ${(c^n)}_{n\in \mathrm{\mathbb{N}}}\subset L^{\infty }(\mathrm{\Omega })$ converging to c in the norm topology of $L^q$ , $q<\infty$ , such that

$$\begin{array}{ccc}\multicolumn{1}{c}{\mathcal{E}[c]\ge \underset{n\to \infty }{\lim \sup }{\mathcal{E}}^n[c^n].}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(37)

Proof: Statement (36) is a direct consequence of the weak lower semi-continuity (32) of the energy functional $\mathcal{E}$ , combined with Lemma 10. Indeed, for any $c^n\rightharpoonup c$ , we have

$$\begin{array}{ccccccc}\multicolumn{1}{c}{\mathcal{E}[c]}& \multicolumn{1}{c}{\le }& \multicolumn{1}{c}{\underset{n\to +\infty }{\lim \inf }\mathcal{E}[c^n]}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{\le }& \multicolumn{1}{c}{\underset{n\to \infty }{\lim \inf }{\mathcal{E}}^n[c^n]+\underset{n\to \infty }{\lim }(\mathcal{E}[c^n]-{\mathcal{E}}^n[c^n])}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{=}& \multicolumn{1}{c}{\underset{n\to \infty }{\lim \inf }{\mathcal{E}}^n[c^n].}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

To prove (37), let us fix $c\in L^{\infty }(\mathrm{\Omega })$ and construct $c^n\in Y_0^n$ as the sequence of piecewise constant approximates of c. By standard approximation theory, we then have $c^n\to c$ in the norm topology of $L^q(\mathrm{\Omega })$ for any $q<\infty$ . Then, the strong continuity of the energy functional $\mathcal{E}$ with respect to the $L^2$ -topology established in Lemma 8 gives

$$\begin{array}{ccc}\multicolumn{1}{c}{\mathcal{E}[c]=\underset{n\to +\infty }{\lim }\mathcal{E}[c^n].}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Moreover, Lemma 9 gives $\mathcal{E}[c^n]\ge {\mathcal{E}}^n[c^n]$ , which directly implies (37). ■

With Theorem 11, it is straightforward to construct global minimizers of the continuum energy (5) as limits of sequences of minimizers of the discrete energy functional (27) on equilateral triangulations.

Theorem 12. Let ${\mathcal{C}}^n\in {\mathrm{ℝ}}^{|{\mathbb{E}}^n|}$ , $n\in \mathrm{\mathbb{N}}$ , be a sequence of global minimizers of the rescaled discrete energy functionals ${\overline{E}}^n$ given by (27). Assume that $\mathrm{su}{p}_{n\in \mathrm{\mathbb{N}}}{\parallel {\mathcal{C}}^n\parallel }_{{\ell }^{\infty }}<\infty$ . Then, up to an eventual selection of a subsequence, $c^n:={\mathbb{Z}}^n[{\mathcal{C}}^n]$ converges weakly-* to $c\in L^{\infty }(\mathrm{\Omega })$ . Under the additional assumption that $c\in C(\overline{\mathrm{\Omega }})$ , c is a global minimizer of the continuum energy functional (5).

Proof: The uniform boundedness of ${\parallel {\mathcal{C}}^n\parallel }_{{\ell }^{\infty }}$ obviously implies that $c^n:={\mathbb{Z}}^n[{\mathcal{C}}^n]$ is uniformly bounded in $L^{\infty }(\mathrm{\Omega })$ and (an eventual subsequence) converges weakly-* to some $c\in L^{\infty }(\mathrm{\Omega })$ . Then, statement (36) of Theorem 11 gives

$$\begin{array}{ccc}\multicolumn{1}{c}{\mathcal{E}[c]\le \underset{n\to +\infty }{\lim \inf }{\mathcal{E}}^n[c^n].}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

We claim that if $c\in C(\overline{\mathrm{\Omega }})$ , then it is a global minimizer of the energy functional (5). For contradiction, let us assume that there exists $\tilde{c}\in C(\overline{\mathrm{\Omega }})$ such that

$$\begin{array}{ccc}\multicolumn{1}{c}{\mathcal{E}[\tilde{c}]\mathcal{E}}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

(38)

By averaging $\tilde{c}$ over the diamonds ${◇}_{\mathrm{ij}}^n$ , we construct a sequence ${({\tilde{\mathcal{C}}}^n)}_{n\in \mathrm{\mathbb{N}}}$ such that ${\tilde{c}}^n:={\mathbb{Q}}^n[{\tilde{\mathcal{C}}}^n]$ converges to $\tilde{c}$ in the norm topology of $L^{\gamma }(\mathrm{\Omega })$ . Moreover, since $\tilde{c}\in C(\overline{\mathrm{\Omega }})$ , the sequence ${\tilde{\mathcal{C}}}^n$ satisfies the property (30). Corollary 7 and the statement (37) of Theorem 11 gives then

$$\begin{array}{ccc}\multicolumn{1}{c}{\underset{n\to \infty }{\lim \sup }{\overline{E}}^n[{\stackrel{\mathrm{~}}{\mathcal{C}}}^n]=\underset{n\to \infty }{\lim \sup }{\mathcal{E}}^n[{\stackrel{\mathrm{~}}{c}}^n]\le \mathcal{E}[\stackrel{\mathrm{~}}{c}].}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

On the contrary, we have

$$\begin{array}{ccc}\multicolumn{1}{c}{\underset{n\to \infty }{\lim \inf }{\overline{E}}^n[{\mathcal{C}}^n]=\underset{n\to \infty }{\lim \inf }{\mathcal{E}}^n[c^n]\ge \mathcal{E}[c].}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

Since by construction we have $\overline{E}[{\mathcal{C}}^n]\le \overline{E}[{\tilde{\mathcal{C}}}^n]$ for all $n\in \mathrm{\mathbb{N}}$ , we finally obtain

$$\begin{array}{ccc}\multicolumn{1}{c}{\mathcal{E}[c]\le \underset{n\to \infty }{\lim \inf }{\overline{E}}^n[{\mathcal{C}}^n]\le \underset{n\to \infty }{\lim \sup }{\overline{E}}^n[{\stackrel{\mathrm{~}}{\mathcal{C}}}^n]\le \mathcal{E}[\stackrel{\mathrm{~}}{c}],}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

which is a contradiction to (38). ■

Finally, let us note that by the standard theory of convex gradient flows on Hilbert spaces, see, e.g., [26,27], we can construct unique weak solutions of the systems (1)–(2). We refer to [1, Section 3.1] for further details.

Proposition 13. Let $r>0$ , $S\in L^2(\mathrm{\Omega })$ satisfying (14) and $c^0\in L^2(\mathrm{\Omega })$ non-negative almost everywhere in Ω. Then, the problems (1) – (4) admit a unique global weak solution $c\in H^1((0,\infty );L^2(\mathrm{\Omega }))$ , with $E_0[c]\in L^{\infty }(0,\infty )$ and satisfying the energy inequality

$$\begin{array}{ccc}\multicolumn{1}{c}{\mathcal{E}[c(t)]+{\int }_0^t{\int }_{\mathrm{\Omega }}{\left|\frac{\mathrm{\partial c}}{\mathrm{\partial t}}(s,x)\right|}^2\mathrm{d}\mathrm{x}\mathrm{d}\mathrm{s}\le \mathcal{E}[c^0].}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

3. Numerical results

In this section, we briefly discuss some of the details of our numerical implementation. We first introduce the finite element spaces, the semi-discrete nonlinear equations, and their time-advancing scheme and then discuss the linearized equations. We conclude with numerical experiments in Section 3.1 and 3.2. For additional details, see [28].

While we used homogeneous Neumann boundary conditions to derive the continuum energy, here we consider the system of equations (1)–(2) closed with more general boundary conditions

$$\begin{array}{ccc}\multicolumn{1}{c}{\{\begin{array}{ccc}\multicolumn{1}{c}{u}& \multicolumn{1}{c}{=g_D\mathrm{on}{\mathrm{\Gamma }}_D,}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{c\frac{\mathrm{\partial u}}{\mathrm{\partial n}}}& \multicolumn{1}{c}{=g_N\mathrm{on}{\mathrm{\Gamma }}_N,}\end{array}}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

where ${\mathrm{\Gamma }}_D\cup {\mathrm{\Gamma }}_N=\partial \mathrm{\Omega }$ , with ${\mathrm{\Gamma }}_D$ (respectively ${\mathrm{\Gamma }}_N$ ) being the subset of the boundary where we may consider non-homogeneous essential (respectively natural) boundary conditions with datum $g_D$ (respectively $g_N$ ). We use finite element spaces consisting of piecewise constant elements for the conductivity and bilinear elements for the pressure in two dimensions, i.e.,

$$\begin{array}{cccc}\multicolumn{1}{c}{C_h}& \multicolumn{1}{c}{=\{c_h\in L^2(\mathrm{\Omega })|{c_h}_{|K}\mathrm{constant,}\mathrm{\forall }K\in \mathrm{\Omega }\},}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{U_h}& \multicolumn{1}{c}{=\{u_h\in H_{{\mathrm{\Gamma }}_D}^1(\mathrm{\Omega })|u_h\in C^0(\mathrm{\Omega }),u_h=g_D\mathrm{on}{\mathrm{\Gamma }}_D,\mathrm{and}{u_h}_{|K}\mathrm{bilinear,}\mathrm{\forall }K\in \mathrm{\Omega }\},}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

where K is a quadrilateral cell of the tessellated domain Ω.

For the sake of stability and well-posedness, we will use a regularized metabolic cost when $\gamma <1$

$$\begin{array}{ccc}\multicolumn{1}{c}{{(c^2+\epsilon )}^{\gamma ∕2},\epsilon \ge 0,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

that reduces to $c^{\gamma }$ when $\epsilon =0$ , as initially proposed in Astuto et al. [29].

After a simple rescaling to ensure the symmetry of the linearized equations, the variational formulation of the discrete problem reads: given $S\in L^2(\mathrm{\Omega })$ , find $u_h(t)\in U_h$ , $c_h(t)\in C_h$ such that, for almost every $t\in (0,T)$

$$\begin{array}{ccccccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\frac{1}{2}{\int }_{\mathrm{\Omega }}\left(\frac{\partial c_h}{\mathrm{\partial t}}{b}_h-|\mathrm{\nabla }{u}_h{|}^2{b}_h+\nu {(c_h^2+\epsilon )}^{(\gamma -2)∕2}{c}_h{b}_h\right)\mathrm{d}\mathrm{x}=0,}& \multicolumn{1}{c}{\mathrm{\forall }{b}_h\in C_h,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{-{\int }_{\mathrm{\Omega }}(c_h+r)\mathrm{\nabla }{u}_h\cdot \mathrm{\nabla }{q}_h\mathrm{d}\mathrm{x}=-{\int }_{\mathrm{\Omega }}S{q}_h\mathrm{d}\mathrm{x}-{\int }_{{\mathrm{\Gamma }}_N}{g}_N{q}_h\mathrm{d}\mathrm{s},}& \multicolumn{1}{c}{\mathrm{\forall }{q}_h\in U_h,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

which leads to the set of differential-algebraic equations

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\frac{1}{2}{\int }_{\mathrm{\Omega }}\left(\frac{\partial c_h}{\mathrm{\partial t}}{\psi }_i-|\mathrm{\nabla }{u}_h{|}^2{\psi }_i+\nu {(c_h^2+\epsilon )}^{(\gamma -2)∕2}{c}_h{\psi }_i\right)\mathrm{d}\mathrm{x}=0,}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{-{\int }_{\mathrm{\Omega }}(c_h+r)\mathrm{\nabla }{u}_h\cdot \mathrm{\nabla }{\phi }_i\mathrm{d}\mathrm{x}=-{\int }_{\mathrm{\Omega }}S{\phi }_i\mathrm{d}\mathrm{x}-{\int }_{{\mathrm{\Gamma }}_N}{g}_N{\phi }_i\mathrm{d}\mathrm{s},}\end{array}$$

where ${\phi }_i$ (respectively ${\psi }_i$ ) are the test functions for $U_h$ (respectively $C_h$ ).

To advance in time, we use the backward Euler scheme, which allows us to preserve the positivity of the conductivity. Dropping the h subscript, we first compute the initial pressure $u_0$ by solving

$${\int }_{\mathrm{\Omega }}(c_0+r)\mathrm{\nabla }{u}_0\cdot \mathrm{\nabla }{\phi }_i\mathrm{d}\mathrm{x}={\int }_{\mathrm{\Omega }}S{\phi }_i\mathrm{d}\mathrm{x}+{\int }_{{\mathrm{\Gamma }}_N}{g}_N{\phi }_i\mathrm{d}\mathrm{s},$$

where $c_0$ is the given initial condition for the conductivity. Denoting by $c_n$ and $u_n$ the finite element approximations for the conductivity and pressure at step n, we then solve the equations

$$\begin{array}{ccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{\frac{1}{2}{\int }_{\mathrm{\Omega }}\left(\frac{c_{n+1}-c_n}{\delta t_n}{\psi }_i-|\mathrm{\nabla }{u}_{n+1}{|}^2{\psi }_i+\nu {(c_{n+1}^2+\epsilon )}^{(\gamma -2)∕2}{c}_{n+1}{\psi }_i\right)\mathrm{d}\mathrm{x}=0,}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{-{\int }_{\mathrm{\Omega }}(c_{n+1}+r)\mathrm{\nabla }{u}_{n+1}\cdot \mathrm{\nabla }{\phi }_i\mathrm{d}\mathrm{x}=-{\int }_{\mathrm{\Omega }}S{\phi }_i\mathrm{d}\mathrm{x}-{\int }_{{\mathrm{\Gamma }}_N}{g}_N{\phi }_i\mathrm{d}\mathrm{s},}\\ \multicolumn{1}{c}{}\end{array}$$

where $\delta t_n$ is the time step.

This nonlinear system of equations is solved using the inexact Newton method [30]. The symmetric indefinite Jacobian matrix at a given linearization point is of the form

$$\begin{array}{ccc}\multicolumn{1}{c}{J=\left[\begin{array}{cc}\multicolumn{1}{c}{A}& \multicolumn{1}{c}{B^T}\\ \multicolumn{1}{c}{B}& \multicolumn{1}{c}{-C}\end{array}\right],}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

where

$$\begin{array}{cccc}\multicolumn{1}{c}{}& \multicolumn{1}{c}{A_{\mathrm{ij}}=\frac{1}{2}{\int }_{\mathrm{\Omega }}(\frac{1}{\delta t_n}+\alpha (c_{n,k})+\beta (c_{n,k})){\psi }_i{\psi }_j\mathrm{d}\mathrm{x},}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{B_{\mathrm{ij}}=-{\int }_{\mathrm{\Omega }}{\psi }_j\left(\mathrm{\nabla }{u}_{n,k}\cdot \mathrm{\nabla }{\phi }_i\right)\mathrm{d}\mathrm{x},}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{}& \multicolumn{1}{c}{C_{\mathrm{ij}}={\int }_{\mathrm{\Omega }}(c_h+r)\mathrm{\nabla }{\phi }_i\cdot \mathrm{\nabla }{\phi }_j\mathrm{d}\mathrm{x},}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

with $c_{n,k}$ (respectively $u_{n,k}$ ) is the linearization point at the kth Newton step for the conductivity (respectively pressure), and the function coefficients α and β are given by

$$\alpha (c)=\nu {(c^2+\epsilon )}^{(\gamma -2)∕2},\beta (c)=\nu (\gamma -2){(c^2+\epsilon )}^{\frac{\gamma -4}{2}}{c}^2.$$

Note that, when $\gamma \ge 1$ and $\epsilon \ge 0$ , A is a symmetric positive definite matrix because

$$\begin{array}{ccc}\multicolumn{1}{c}{\frac{1}{\delta t_n}+\alpha (c)+\beta (c)\ge \nu {(c^2+\epsilon )}^{(\gamma -4)∕2}(\epsilon +(\gamma -1)c^2)\ge 0.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

When $\gamma <1$ , the positive definiteness of A cannot be guaranteed unless we make assumptions on the time step.

The matrix J admits the exact factorization

$$\left[\begin{array}{cc}\multicolumn{1}{c}{A}& \multicolumn{1}{c}{B^T}\\ \multicolumn{1}{c}{B}& \multicolumn{1}{c}{-C}\end{array}\right]=\left[\begin{array}{cc}\multicolumn{1}{c}{A}& \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{B}& \multicolumn{1}{c}{I}\end{array}\right]\left[\begin{array}{cc}\multicolumn{1}{c}{A^{-1}}& \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}& \multicolumn{1}{c}{-D}\end{array}\right]\left[\begin{array}{cc}\multicolumn{1}{c}{A}& \multicolumn{1}{c}{B^T}\\ \multicolumn{1}{c}{0}& \multicolumn{1}{c}{I}\end{array}\right],D=C+B{A}^{-1}{B}^T,$$

(39)

where the Schur complement D can be explicitly assembled since A is a block-diagonal matrix with each block associated with a single element. The inverse factorization requires solving the inexpensive diagonal problem A twice, and the Schur complement system D only once. Numerically, we observed that D is a symmetric positive semi-definite matrix representing a perturbed scalar diffusion problem having the same kernel as the discretized Poisson problem C. Instead of solving the Schur complement system exactly, we can thus use an algebraic multigrid (AMG) preconditioner, which guarantees robustness and has an overall cost that is linear in the matrix size.

The code used to conduct the numerical experiments described in this section is publicly available as a tutorial of the Portable and Extensible Toolkit for Scientific Computing (PETSc, [31]). ¹ The management of unstructured grids and the implementation of finite elements are based on the DMPLEX infrastructure [32]. Time integration is performed using the TS module [33]; adaptive time stepping is performed using digital filtering-based methods [34] combined with the approximation of the local truncation error using extrapolation. The nonlinear systems of equations arising from the time-advancing schemes are solved with a tight absolute tolerance of 1.E-10; Jacobian linear systems are solved using the right-preconditioned generalized minimal residual (GMRES) method [35], and the Schur complement system is approximatively solved using the application of an AMG preconditioner based on smoothed aggregation [36]. Linear system tolerances are dynamically adjusted using the well-known Eisenstat and Walker trick, which prevents over-solving in the early steps of the Newton process and tightens the accuracy as convergence is approached [37]. The nonlinear time step problems can sometimes be quite challenging to solve; as a safeguard, we cap the maximum number of nonlinear iterations at 10 and shorten the time step if the Newton process stagnates and does not reach convergence within the maximum number of iterations.

3.1. Network formation

In the first numerical experiment, we demonstrate that the scalar model is capable of simulating a network formation process. For this experiment, we consider a unit square domain $\mathrm{\Omega }={[0,1]}^2$ , homogeneous Neumann boundary conditions (i.e. ${\mathrm{\Gamma }}_D=\varnothing$ , $g_N=0$ ), and the model parameters $r=0.0001$ , $\epsilon =0.001$ , $\gamma =0.75$ , and $\nu =0.05$ . The domain is discretized using a structured quadrilateral mesh with a 512 x 512 cell arrangement, resulting in a total of 263,000 degrees of freedom. The initial conductivity field is the constant function $c_0=1$ , the simulation is run up to the final time $T=200$ , and we use the Gaussian source

$$S(x)=S_0(x)-{\int }_{\mathrm{\Omega }}{S}_0(x)\mathrm{d}\mathrm{x},S_0(x)=e^{-500{\parallel x-x_0\parallel }^2},$$

with $x_0=(0.25,0.25)$ . Figure 1 presents the system energy E (left panel), the time step size $\delta t_n$ (center), and the number of nonlinear iterations (right) as functions of the simulated time. The legend in the right panel reports the total number of time steps, the total number of Newton steps, and the average number of Krylov iterations per Newton step, providing an overview of the algorithmic costs of the solver. Snapshots of the conductivity field $c_h$ (in logarithmic scale) taken at selected times indicated by the text boxes in Figure 1, are shown in Figure 2.

Figure 1.

Network formation: energy (left panel, logarithmic scale), time step (central, logarithmic scale), and number of nonlinear iterations as a function of simulation time for 512 x 512 mesh.

Figure 2.

Network formation: conductivity in logarithmic scale at selected time instances as reported at the top of each panel (see Figure 1).

The system energy decreases monotonically throughout the simulation. In the initial phase, the energy decays rapidly as the conductivity field evolves from an initially homogeneous state to a spatially heterogeneous configuration (see panel (a) of Figure 2), allowing relatively large time steps. At this stage, the dynamics are dominated by the minimization of the term $(c+r)|\mathrm{\nabla }u{|}^2$ . In the subsequent phase, network formation begins (see panels (b) and (c)), driven by the interplay between this positive term and the negative contribution from the metabolic term. During this process, energy variations become smaller, and the solver selects smaller time steps. In the final phase, the energy decreases more slowly toward its steady state, with minimal changes in the conductivity field (see panels (d) and (e)). Overall, the problem exhibits a mild nonlinearity. The backward Euler time advancing scheme encounters a few nonlinear solver failures during the later stages of network formation. Nevertheless, the average number of Krylov iterations per Newton step remains very low, confirming the effectiveness and robustness of the linear preconditioning strategy.

3.2. Solutions of the p-Laplacian equation

We dedicate a substantial portion of the numerical section to the integration of the differential-algebraic equations (1) and (2) up to their steady state, in order to compute solutions of the p-Laplacian equation

$$\begin{array}{ccc}\multicolumn{1}{c}{-\mathrm{\nabla }\cdot (|\mathrm{\nabla }u{|}^{p-2}\mathrm{\nabla }u)=S,p=\frac{2\gamma }{\gamma -1}>2,}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

with $\gamma >1$ , $r=0$ , $\epsilon =0$ , and $\nu =1$ . To evaluate the approximation properties of the method, we first apply the method of manufactured solutions and compare the numerical results against analytical solutions. Finally, we present experiments demonstrating the robustness of the method for large values of p and non-convex domains.

A priori estimates for the lowest order $H^1$ -conforming finite element discretizations of the p-Laplacian equation with $p>2$ have been established in Ebmeyer and Liu [38] (see also Diening and Ruzicka [39]) in terms of a quasi-norm for convex domains and a sufficiently smooth source term S

$$|u-u_h{|}_{(u,p)}^2\le C{h}^2{\int }_{\mathrm{\Omega }}|\mathrm{\nabla }u{|}^{p-2}|D^{2}u{|}^2\mathrm{d}\mathrm{x},$$

(40)

where the quasi-norm $|\cdot {|}_{(w,p)}$ is defined as

$$|v{|}_{(w,p)}^2={\int }_{\mathrm{\Omega }}(|\mathrm{\nabla }w|+|\mathrm{\nabla }v|{)}^{p-2}|\mathrm{\nabla }v{|}^2\mathrm{d}\mathrm{x}.$$

(41)

In addition to the quasi-norm, we report convergence rates in terms of the standard Lebesgue $L^p(\mathrm{\Omega })$ norm and the semi-norm of the Sobolev space $W^{1,p}(\mathrm{\Omega })$

$$\begin{array}{ccc}\multicolumn{1}{c}{\parallel v{\parallel }_p=({\int }_{\mathrm{\Omega }}|v{|}^p\mathrm{d}\mathrm{x}{)}^{1∕p},|v{|}_p=({\int }_{\mathrm{\Omega }}|\mathrm{\nabla }v{|}^p\mathrm{d}\mathrm{x}{)}^{1∕p}.}& \multicolumn{1}{c}{}& \multicolumn{1}{c}{}\end{array}$$

We consider three test cases from Barrett and Liu [40] with radially symmetric smooth functions on the unit square. The first two tests use

$$\begin{array}{cc}\begin{array}{cc}\multicolumn{1}{c}{u(r)}& \multicolumn{1}{c}{=\frac{p-1}{\sigma +p}{(\sigma +2)}^{\frac{1}{1-p}}(1-r^{\frac{\sigma +p}{p-1}}),}\\ \multicolumn{1}{c}{S(r)}& \multicolumn{1}{c}{=r^{\sigma },}\\ \multicolumn{1}{c}{}\end{array}& \end{array}$$

(42)

with $\sigma =0,p=4$ (denoted by TC1) and $\sigma =7,p=4$ (TC2) as given in Examples 5.2 and 5.3 in Barrett and Liu [40]. The third test case (denoted by TC3) is also from Barrett and Liu [40] (Example 5.4)

$$\begin{array}{c}\begin{array}{cc}\multicolumn{1}{c}{u(r)}& \multicolumn{1}{c}{=\{\begin{array}{cc}\multicolumn{1}{c}{0}& \multicolumn{1}{c}{r<a}\\ \multicolumn{1}{c}{{(r-a)}^4}& \multicolumn{1}{c}{r\ge a}\end{array}}\\ \multicolumn{1}{c}{S(r)}& \multicolumn{1}{c}{=\{\begin{array}{cc}\multicolumn{1}{c}{0}& \multicolumn{1}{c}{r<a}\\ \multicolumn{1}{c}{4^{p-1}{(r-a)}^{3p-4}(2-3p+\frac{a}{r})}& \multicolumn{1}{c}{r\ge a}\end{array}}\end{array}\end{array}$$

(43)

with $a=0.3$ and $p=4$ . In order to showcase the convergence of the method for larger p-Laplacian exponents, we complement these tests with a fourth one, denoted by TC4, considering (42) with $\sigma =7$ and $p=20$ . For all these tests, we use non-homogeneous essential boundary conditions defined by the analytical solution at the boundary, i.e., ${\mathrm{\Gamma }}_D=\mathrm{\Gamma }$ and ${\mathrm{\Gamma }}_N=\varnothing$ .

We consider five levels of uniform refinement starting from a uniform grid of $16\times 16$ elements and report the convergence curves of the final pressures in Figure 3. Optimal first-order convergence in the quasi-norm $|\cdot {|}_{(u,p)}$ is obtained in all the cases. Convergence rates in $L^p(\mathrm{\Omega })$ and $W^{1,p}(\mathrm{\Omega })$ instead depend on the test case. For TC1, we recover the same convergence rates as in Barrett and Liu [40] (see Table 5.2 in the reference), while our convergence rates for TC2 and TC3 are different than those presented in Tables 5.3 and 5.4 in Barrett and Liu [40] because Barrett and Liu reported the error $\parallel {\pi }_{h}u-u_h\parallel$ , where ${\pi }_{h}u$ is the finite element method (FEM) interpolant of the exact solution, thus observing superconvergence. The convergence rates in $L^p(\mathrm{\Omega })$ and $W^{1,p}(\mathrm{\Omega })$ reported in Figure 3 are in perfect agreement with the best approximation error for linear quadrilateral elements (see, e.g., Corollary 4 in Ainsworth and Kay [41]). For TC1, we have that $u\in W^{m,p}(\mathrm{\Omega })$ with $m=\frac{4}{3}+\frac{1}{2}=\frac{16}{11}\approx 1.83$ , while for TC4 $m=\frac{27}{19}+\frac{2}{20}=\frac{289}{190}\approx 1.52$ . For TC2 and TC3, we have $m>2$ , and thus, we recover optimal convergence rates using bi-linear elements.

Figure 3.

MMS errors and convergence rates (in parentheses) for the test cases TC1, TC2, TC3, and TC4 and different error metrics.

Details on the solution algorithm for the cases TC2 and TC4 are reported in Figure 4; for each time step n, we report, from left to right, the time increment $\delta t_n$ , the p-Laplacian energy $\frac{1}{p}{\int }_{\mathrm{\Omega }}|\mathrm{\nabla }{u}_h{|}^p-S{u}_h\mathrm{d}\mathrm{x}$ , the $L^2$ norm of $|\mathrm{\nabla }{u}_n{|}^2-c_n^{\gamma -1}$ , and the number of Newton steps as a function of simulation time for the different levels of refinement, from the coarsest $r_0$ , to the finest $r_4$ . The legends in the right panels display the total number of time steps, the total number of Newton steps, and the average number of Krylov iterations per Newton step, providing additional information on the computational costs of the method.

Figure 4.

From left to right: time step, p-Laplacian energy $1∕p{\int }_{\mathrm{\Omega }}|\mathrm{\nabla }{u}_h{|}^p-S{u}_h$ (final value in the legend), $L^2$ norm of $|\mathrm{\nabla }{u}_n{|}^2-c_n^{\gamma -1}$ , and number of Newton steps as a function of time for different levels of uniform refinement r for TC2 (top row) and TC4 (bottom row).

In both test cases, the time step is adaptively adjusted to follow the dynamics of the differential-algebraic equations, which is smoother in the TC2 case. The p-Laplacian energy and $\parallel |\mathrm{\nabla }{u}_n{|}^2-c_n^{\gamma -1}{\parallel }_2$ decrease over time, with $\parallel |\mathrm{\nabla }{u}_n{|}^2-c_n^{\gamma -1}{\parallel }_2$ being reduced up to the discretization error $\mathcal{O}(h)$ because we evaluate the quantity at quadrature points and $\mathrm{\nabla }{u}_n$ is not piecewise constant when using bi-linear quadrilateral elements for $U_h$ . For TC2, the number of time steps and the total number of Newton steps remain independent of the refinement level, making the convergence of the method mesh-independent. The linear preconditioning strategy is also highly effective, requiring, on average, only two Krylov iterations per Newton step in all cases. For the TC4 case, time relaxation is more challenging, and a larger value of p leads to stiffer dynamics at all refinement levels. After an initial transient phase, during which the Poisson problem with constant coefficients relaxes smoothly, the increasing nonlinearity of the problem causes the method to enter a second phase that demands a higher number of both time and Newton steps. This occurs because we cannot guarantee the non-negativity of $c_h$ within the nonlinear time step solver, which in turn fails more frequently, resulting in the time step being forced to shrink. Nevertheless, as the mesh is refined, these numbers appear to plateau. To address these challenges, we are currently investigating extensions of our method based on a mass-preserving mixed formulation of the Poisson problem [42]. We are also exploring alternative time-advancing schemes such as pseudo-transient continuation [43], where the rule used to adapt the time step differs fundamentally from the conventional local error–control strategies employed in ODE solvers. This distinction is natural, as the objective of reaching equilibrium as efficiently as possible contrasts with that of accurately tracking a specific trajectory over time.

We now consider a more challenging test case with a non-convex domain and mixed boundary conditions for the pressure, using an analytical solution expressed in polar coordinates

$$\begin{array}{cc}\begin{array}{cc}\multicolumn{1}{c}{u(r,\theta )}& \multicolumn{1}{c}{=r^{\alpha }\sin (\mathrm{\alpha \theta }),}\\ \multicolumn{1}{c}{S(r,\theta )}& \multicolumn{1}{c}{=-\alpha |\alpha {|}^{p-2}(\alpha -1)(p-2)r^{\mathrm{p\alpha }-\alpha -p}\sin (\mathrm{\alpha \theta }),}\\ \multicolumn{1}{c}{}\end{array}& \end{array}$$

(44)

with $\alpha =2$ and $p=4$ . This test case, referred to as TC5, is defined on an L-shaped domain $\mathrm{\Omega }={(-1,1)}^2\setminus [0,1)\times (-1,0]$ . We impose non-homogeneous essential boundary conditions for the pressure on ${\mathrm{\Gamma }}_D=(\{0\}\times [-1,0])\cup ([0,1]\times \{0\})$ where the prescribed values are taken from the analytical solution. On the complementary part of the boundary ${\mathrm{\Gamma }}_N=\mathrm{\Gamma }\setminus {\mathrm{\Gamma }}_D$ , we enforce non-homogeneous natural boundary conditions with $g_N=|\mathrm{\nabla }u{|}^{p-2}\frac{\mathrm{\partial u}}{\mathrm{\partial n}}$ where u is the exact solution. Convergence rates with different levels of uniform refinement are given in Figure 5, showing second-order convergence in the $L^p$ norm and superlinear convergence in the $W^{1,p}(\mathrm{\Omega })$ semi-norm and in the quasi-norm $|\cdot {|}_{(u,p)}$ .

Figure 5.

TC5: MMS errors and convergence rates (in parentheses) for different error metrics.

We conclude this section by presenting results for different values of p in a challenging problem without an analytical solution, following the experimental setting of Section 7.3 in Balci et al. [23]. This test, referred to as TC6, considers a constant source term $S=2$ and homogeneous essential boundary conditions for the pressure on the L-shaped domain ${(-1,1)}^2\setminus [0,1)\times (-1,0]$ . The domain is uniformly discretized using 10,000 elements, with mesh size $h=1.9\times 10^{-3}$ . Solver performance for p-Laplacian exponents $p=5,10,20,50$ (corresponding to $\gamma =5∕3,5∕4,10∕9,25∕24$ ) is reported in Figure 6.

Figure 6.

TC6: time step (left panel), $L^2$ norm of $|\mathrm{\nabla }{u}_n{|}^2-c_n^{\gamma -1}$ (center), and the number of Newton steps (right) as a function of time for different p-Laplacian exponents p.

The final values of $\parallel |\mathrm{\nabla }{u}_n{|}^2-c_n^{\gamma -1}{\parallel }_2$ are larger for larger p, possibly because we did not use a p-dependent norm. The solver is not p-independent overall, as it requires more time and Newton steps as p is increased, due to the increase in nonlinear solver failures. However, the average number of Krylov iterations per Newton step remains essentially constant across all cases, confirming the robustness of the Schur complement-based preconditioner. Figure 7 shows the final states of the pressure $u_h$ and the conductivity $c_h$ for $p=50$ , together with the pointwise error $|\mathrm{\nabla }{u}_h{|}^2-c_h^{\gamma -1}$ (evaluated in post-processing) and its projection onto $L^2(\mathrm{\Omega })$ , ${\pi }_h(|\mathrm{\nabla }{u}_h{|}^2-c_h^{\gamma -1})$ (evaluated within the finite element simulation code). The conductivity (plotted on a logarithmic scale) is small in regions of sharp pressure gradients, where the singularities appear to have Hausdorff dimension one, and attains large values near the concave corner of the domain. These patterns closely resemble the mesh adaptation patterns reported in Figure 7.3 of Balci et al. [23]. The error distribution is nearly uniform across the domain and noticeably larger near the L-shaped corner; the $L^2$ projection of the error is small everywhere in the domain and slightly larger near the corner and where $c_h\approx 0$ . We remark that a small regularization parameter, $\epsilon =10^{-8}$ , was introduced to ensure convergence of the solver for $p=20$ and $p=50$ . This small level of regularization does not affect the numerical approximation properties, as the region of the domain where $c_h<\epsilon$ has Hausdorff dimension one. Future work will focus on adaptive mesh refinement strategies and the development of error indicators tailored explicitly to the differential-algebraic transient dynamics of the solver.

Figure 7.

Top row: final state for pressure (left) and conductivity (right) for the TC6 test case with $p=50$ (corresponding to $\gamma =25∕24$ ).