1 Introduction↩︎

This work concerns the numerical solution of minimization problems of the type \[\label{eq:energy95general} \mathop{\mathrm{argmin}}_{u\in W}\int_{\Omega} \phi(|Du|)\,\mathrm{d}x- F(u),\tag{1}\] where \(\Omega\) is an open, bounded domain \(\Omega\subset\mathbb{R}^d\), \(d=2,3\), the integrand is governed by an N-function \(\phi\) [1], and \(D\) denotes either the gradient of a scalar field or the symmetric gradient \(\varepsilon({\boldsymbol{u}}) = \frac{1}{2}(\nabla{\boldsymbol{u}}+ \nabla{\boldsymbol{u}}^T)\) of an incompressible vector field in an Orlicz-Sobolev functional space \(W\). Such energies arise in nonlinear diffusion, optimal design, and in the modeling of generalized Newtonian fluids. The associated Euler–Lagrange equations are nonlinear and may be degenerate, which presents challenges for both analysis and numerical approximation.

The finite element approximation of problems of this type has a long history, beginning with classical works on the \(p\)-Laplacian and on generalized Newtonian models such as Carreau and power-law (aka Ostwald and de Waele fluids) [2]–[6]. FEM optimality was shown in [7] in terms of a quasi-norm; adaptive finite element methods for the \(p\)-Laplacian were analyzed in [8]–[10], while the finite element approximation of \(p\)-Stokes systems was studied in [11]. High-order FEM approximation properties for the \(p\)-Laplacian can be found in [12]. Discretizations of energies with non-standard growth have also been investigated, for instance in [13]. Structurally related developments include discontinuous Galerkin methods for quasilinear elliptic problems [14], [15], and, more recently, Hybrid High-Order methods [16]–[18] and Virtual Element methods [19], [20].

The nonlinear algebraic systems arising from these discretizations are commonly solved either by Newton-type methods or by fixed-point linearizations. Among the latter, the Kačanov method is one of the most widely used approaches. It originated in nonlinear magnetostatics [21] and was later developed as a general linearization technique for nonlinear PDEs and variational inequalities [22]–[24], see also [25]. For energies of the type 1 , a relaxed version of the Kačanov method has been studied for the \(p\)-Laplacian [26], [27], shear-thinning Carreau-type fluids [28], and regularized Bingham models [29]. For a damped version of the algorithm, see [30], and for recent work on the relaxed version, see [31]. Finally, closely related to the Kačanov method is the adaptive iterative linearized Galerkin method [32].

The starting point of the present work is the differential-algebraic system introduced in [33] for modeling evolution of biological transportation networks \[\frac{ds}{dt} = |\nabla u|^2 - s^{\gamma-1}, \quad -\nabla\cdot\bigl((s+r)\nabla u\bigr) = f,\quad \gamma > 1,\] where \(r>0\) is an arbitrarily small regularization parameter ensuring the well-posedness of the linear elliptic problem. The steady states of this system solve the Euler–Lagrange equations of the regularized \(p\)-Laplacian \[\int_\Omega (|\nabla u|^{p-2} + r)\nabla u \cdot \nabla v \,\mathrm{d}x = \int_\Omega f v \,\mathrm{d}x, \quad p = \frac{2\gamma}{\gamma-1}.\] Rewriting \(\gamma\) in terms of \(p\) and introducing the change of variables \(c=s^{2/(p-2)}\), the dynamics become \[\frac{dc}{dt} = \frac{2}{p-2} c^{\frac{4-p}{2}}\bigl(|\nabla u|^2-c\bigr), \quad -\nabla\cdot\left( \left(c^{\frac{p-2}{2}}+r\right) \nabla u\right) = f .\] If the goal is to recover only the steady state, the prefactor in the ordinary differential equation (ODE) is not essential, since it changes the parametrization of the relaxation in time, but not its equilibrium. One is thus tempted to deliberately drop this factor and use the simpler evolution \[\label{eq:ode95intro} \frac{dc}{dt} = |\nabla u|^2 - c.\tag{2}\] This seemingly minor simplification was the starting point of the present paper: it produces an auxiliary dynamics that is no longer tied to the specific algebraic structure of the power-law model (which is only encoded through the elliptic constraint) and can therefore be extended to more general constitutive laws. Note that, for \(r=0\), a forward Euler discretization of 2 with time step equal to one gives the Kačanov iteration for the \(p\)-Laplacian [26], \[c_{n+1} = |\nabla u_n|^2, \quad \int_\Omega c_{n+1}^{(p-2)/2}\,\nabla u_{n+1}\cdot \nabla v \,\mathrm{d}x = \int_\Omega f v \,\mathrm{d}x.\]

The aim of this work is not to derive an equivalent minimization problem for 1 . Rather, we introduce an auxiliary-field formulation in which the nonlinearity is shifted from the state equation to a scalar variable \(c\) representing, at equilibrium, the squared magnitude \(|Du|^2\). The resulting formulation provides a variationally consistent framework for designing and analyzing approximation schemes based on a sequence of uniformly elliptic linear state problems of the type \[\label{eq:constraint95gen} \int_{\Omega}\mu(c(x))\,D u\!:\!D v\,\mathrm{d}x = F(v) ,\quad\mu(r):=\frac{\phi'(\sqrt{r})}{\sqrt{r}}.\tag{3}\] Note that, in the non-Newtonian fluids literature (see e.g. [34], the viscosity is written as a function \(\widehat\mu(\dot{\gamma})\), where the shear rate \(\dot{\gamma}\) is the magnitude of the rate-of-strain tensor \[\dot{\gamma} = \sqrt{2}|\varepsilon({\boldsymbol{u}})|.\] Here the auxiliary variable tracks, modulo the \(\sqrt{2}\) scaling, the squared modulus \(|\varepsilon({\boldsymbol{u}})|^2\).

In the first part of the paper, we construct an auxiliary energy with domain of arbitrarily truncated functions, i.e. \(\epsilon \le c(x) \le 1/\epsilon\) almost everywhere for an arbitrary \(\epsilon \in (0,1)\), and show that, under suitable assumptions, it is lower semicontinuous on a complete geodesic metric space whose structure depends on the first derivative of the viscosity. A central point of the analysis is that the space to which the auxiliary variable \(c\) belongs is isometric to \({L^2(\Omega)}\). We then use the theory of gradient flows in metric spaces [35] to prove exponential convergence, independent of \(\epsilon\), of the method of minimizing movements [36] toward the minimizers in the truncated space. The convergence mechanism is the geodesical \(\lambda\)-convexity of the auxiliary energy in this adapted metric. In particular, for the power-law case, this yields a rigorous convergence theory for \(4/3\le p\le 4\), \(p\ne2\), with \(p=2\) excluded, corresponding to the linear problem.

The second part of the paper connects the continuous formulation with finite element discretizations. Instead of discretizing the method of minimizing movements directly, we derive a finite-dimensional evolution leveraging the theory of metric gradient flows [37] (see also [38]). The discrete metric induces an explicit Riesz map, which is crucial for identifying first variations with gradients, and yields an evolution of the form 2 ; the dependence on the model rheology enters only through the constraint 3 . We prove global existence and uniqueness of the solution of this constrained ODE on the positive cone of the discrete auxiliary space under the weaker assumption that the \(\mu'(r)\) has a fixed sign and is nonzero for \(r>0\). Under additional structural assumptions that match the continuous theory, we also prove convergence to the finite element solution of the Euler–Lagrange equations associated with 1 . In particular, for power-law models this gives a rigorous convergence result in the range \(4/3\le p\le 4\), \(p\ne2\) that matches the continuous analysis. The cases \(p<4/3\) and \(p>4\) are not covered by the present theory, but the same framework suggests asymptotic convergence-rate estimates near equilibrium that are consistent with our numerical experiments.

We then describe several time-discretization strategies for the auxiliary metric gradient flow, which contain the Kačanov iteration as a special case. Motivated by pseudo-transient continuation methods [39] and their extensions to differential-algebraic equations [40], we also propose a time-adaptive pseudo-transient iteration that can be implemented using scalable linear solvers. A convergence analysis of these fully discrete schemes is not pursued here and is left for future work. Nevertheless, the numerical experiments demonstrate the robustness of the methods across different models, confirm the predicted rates for power-law models and for the Carreau–Yasuda model [41], [42] in the theoretically covered range, and support the asymptotic estimates outside that range. Numerical results also show that the methods are mesh-independent and indicate that the proposed strategy is competitive with, and in several cases more robust than, Newton’s method, including solving the \(p\)-Laplacian equations with \(p=100\).

The paper is organized as follows. Section 2 introduces the Orlicz–Sobolev setting and the metric structure used for the auxiliary variable. Section 3 constructs the auxiliary energy in the scalar and incompressible vector-valued cases, and proves lower semicontinuity and geodesical \(\lambda\)-convexity estimates. Section 4 presents the finite element discretization, derives the discrete metric gradient flow, proves well-posedness and convergence of the semi-discrete ODE under suitable assumptions, and introduces the time discretizations. Section 5 reports numerical experiments for power-law, Carreau–Yasuda, regularized Bingham, and optimal-design models. Section 6 concludes with a discussion of open directions.

2 Preliminaries↩︎

2.1 Orlicz-Sobolev spaces↩︎

The function \(\phi'_+\) is non-decreasing, satisfies \(\phi'_+(0) = 0\) and \(\phi'_+(r) > 0\) for \(r>0\), with \(\lim_{r \to +\infty}\phi_+'(r)=+\infty\), and is strictly increasing on \((0,+\infty)\) [1]. For simplicity, we assume throughout that \(\phi\) has enough regularity for the following derivatives to be well defined:

We also assume the following growth condition:

This assumption ensures that \(\phi\) and its complementary function \(\phi^*\) defined as \[\label{eq:phi95complement} \phi^\ast(r):=\int^r_0 (\phi')^{-1}(s)\,\mathrm{d}s,\tag{4}\] satisfy a \(\Delta_2\) condition [44], i.e. \(\exists \,C \geq 1 : \phi(2r)\le C\phi(r),\) for all \(r\ge 0\). Under Assumption 2, the natural spaces to study energies of the form 1 are the Orlicz-Sobolev spaces [45]. Associated with \(\phi\) is the Orlicz space \({L^\phi(\Omega)}\), which is a separable and reflexive Banach space equipped with the Luxemburg norm \[\|v\|_{{L^\phi(\Omega)}} :=\inf\Big\{\lambda>0:\int_{\Omega}\phi(|v|/\lambda)\,\mathrm{d}x\le 1\Big\},\] and, under the sole \(\Delta_2\) condition, we have the equivalence (see e.g. [46]) \[v\in {L^\phi(\Omega)}\iff \int_\Omega \phi(|v|) \,\mathrm{d}x< +\infty.\] The Orlicz–Sobolev space \(W^{1,\phi}(\Omega)\) consists of all \(v\in {L^\phi(\Omega)}\) such that \(\nabla v\in L^\phi(\Omega;\mathbb{R}^d)\), and the following continuous embeddings hold for bounded domains [47]) \[\begin{align} \phi(r)/r^2 \text{ non-decreasing}&\Rightarrow {L^\phi(\Omega)}\hookrightarrow {L^2(\Omega)},\quad W^{1,\phi}(\Omega) \hookrightarrow H^1(\Omega),\tag{5}\\ \phi(r)/r^2 \text{ non-increasing}&\Rightarrow {L^2(\Omega)}\hookrightarrow {L^\phi(\Omega)},\quad H^1(\Omega) \hookrightarrow W^{1,\phi}(\Omega).\tag{6} \end{align}\]

Instead of working with the energy functional 1 , we consider a reparameterization in terms of \(|D u|^2\), and introduce the function \[\label{eq:Phi} \Phi(r):=\phi(\sqrt{r}),\tag{7}\] so that \(\phi(|D u|)=\Phi(|D u|^2)\); from ?? we have the following characterization \[\label{eq:Phi95mu} \Phi(r) := \frac{1}{2}\int^r_0\mu(s)\,\mathrm{d}s,\quad \mu(r):=2\Phi'(r)=\frac{\phi'(\sqrt{r})}{\sqrt{r}}.\tag{8}\] The function \(\Phi\) is continuous, satisfies \(\Phi(0)=0\), and \(\Phi(r)>0\) for \(r>0\) and it is strictly increasing, thus \(\mu(r) > 0\) for \(r > 0\); however, since the mapping \(r\mapsto \sqrt{r}\) is concave, \(\Phi\) is generally not an N-function. A direct calculation yields \[\label{eq:Phipp} \Phi''(r) = \frac{1}{2} \mu'(r) = \frac{1}{4 r\sqrt{r}}\left(\phi''(\sqrt{r})\sqrt{r} - \phi'(\sqrt{r})\right) = \frac{1}{4 t^3}\left(\phi''(t)t - \phi'(t)\right) =\frac{1}{4 t} \frac{d}{dt}\left(\frac{\phi'(t)}{t}\right), \quad t=\sqrt{r}.\tag{9}\] This leads to the next assumption, which we require to hold throughout this work.

From 9 we thus have that \(\Phi\) is strictly convex (resp. concave) if and only if \(\phi'(r)/r\) is increasing (resp. decreasing). and we derive the convenient representation \[\label{eq:mup95exact95final} \mu'(r)=\frac{1}{2r}\left(\frac{t\phi''(t)}{\phi'(t)}-1\right)\mu(r). \quad t = \sqrt{r}.\tag{10}\] Note that the assumed convexity or concavity of \(\Phi\) guarantees, by a simple application of the monotone form of de L’Hôpital’s rule, that one of the embeddings 5 –6 holds. This result is stated without proof in the next Lemma.

For our analysis, we will need two further assumptions. The first is the so-called uniform convexity assumption introduced in [48].

Note that the following chain of inequalities holds \[1 < 1 + {{S}^{-}_{\phi}}\le {{R}^{-}_{\phi}}\le {{R}^{+}_{\phi}}\le 1 + {{S}^{+}_{\phi}}< +\infty.\] We state here some useful Lemmas, for which we provide proofs in the Appendix. The first Lemma provides bounds for two-point growth estimates of \(\phi'\) and \(\mu\); for the proof, see Lemma 7.

Using Lemma 2, we can prove the following; the proof is postponed in Lemma 8.

2.2 Metric space↩︎

Associated with \(\Phi\) is the function space \({L^\Phi(\Omega)}\), defined, with some abuse of notation, as \[\label{eq:LPhi95def} {L^\Phi(\Omega)}:= \left\{ c : \Omega \rightarrow \mathbb{R}\text{ measurable} : \sqrt{|c|} \in {L^\phi(\Omega)}\right\}.\tag{11}\] Note that \({L^\Phi(\Omega)}\) is a reflexive and separable Orlicz space if and only if \(\Phi\) is convex. When \(\Phi\) is concave, the Luxemburg norm is not well defined, and \({L^\Phi(\Omega)}\) is not an Orlicz space. Nevertheless, the modular is always well defined, and we have the equivalences \[\label{eq:c95modular} c\in {L^\Phi(\Omega)} \Longleftrightarrow \int_\Omega \phi \left(\sqrt{|c|}\right)\,\mathrm{d}x<+\infty \Longleftrightarrow \int_\Omega \Phi (|c|)\,\mathrm{d}x<+\infty.\tag{12}\] To prove that \({L^\Phi(\Omega)}\) is a function space also in the concave case, we only need to show it is closed under addition, since the scaling by constants clearly yields functions in \({L^\Phi(\Omega)}\). Since \(\Phi\) is non-decreasing, \(\Phi(0)=0\), and it is subadditive on \([0,+\infty)\), we have \(\Phi(|r_1+r_2|)\le \Phi(|r_1|+ |r_2|)\le \Phi(|r_1|)+\Phi(|r_2|)\), for all \(r_1,r_2\), which, when integrated, gives \(c_1+c_2\in {L^\Phi(\Omega)}\) for all \(c_1,c_2\in{L^\Phi(\Omega)}\).

In order to define a topology on \({L^\Phi(\Omega)}\) which is valid in both the convex and concave case, we introduce a further assumption.

We then introduce the odd function \[\label{eq:psi95def} \psi(r) := \text{sgn}(r)\int^{|r|}_0 \sqrt{\frac{\omega}{2}\mu'(s)}\,\mathrm{d}s,\tag{13}\] where \[\label{eq:omega95def} \omega := \begin{cases} \phantom{-}1, & \Phi\text{ convex},\\ -1, &\Phi\text{ concave}, \end{cases}\tag{14}\] and note that \(\psi\) is a continuous, strictly increasing function under Assumption 3. The following Lemma proves that under Assumption 5 it is a homeomorphism and thus its inverse \(\psi^{-1}(r)\) is also a continuous, strictly increasing function; this allows us to identify functions of \({L^\Phi(\Omega)}\) with functions of \({L^2(\Omega)}\) via the bijective map \[\label{eq:psi95map95def} \Psi : {L^\Phi(\Omega)}\rightarrow {L^2(\Omega)}, \quad \Psi(c)(x) := \psi(c(x)),\tag{15}\] and obtain the characterization \(c \in {L^\Phi(\Omega)}\iff \Psi(c) \in {L^2(\Omega)}\). For the proof, see Lemma 9.

Using the bijective map \(\Psi\) given in 15 we can define a metric on \({L^\Phi(\Omega)}\) as \[\begin{align} \label{eq:X95distance} d(c_0,c_1) := \|\Psi(c_1) - \Psi(c_0)\|_{{L^2(\Omega)}}. \end{align}\tag{16}\] Hence, \(d(\cdot,\cdot)\) is the usual distance in \({L^2(\Omega)}\) under the change of variables \(\eta = \Psi(c)\). Clearly, non-negativity, symmetry, and triangle inequality hold; since \(\Psi\) is a bijection, we can show that \(d(u,v)=0\) implies \(u=v\), and thus \(d(\cdot,\cdot)\) is a metric. A key tool in our continuous analysis is the existence of a geodesic, denoted by \(\gamma_g(s)\), connecting two arbitrary points \(c_0\) and \(c_1\) of \({L^\Phi(\Omega)}\), i.e. \[\gamma_g : [0,1] \rightarrow{L^\Phi(\Omega)},\quad \gamma_g(0) = c_0, \quad \gamma_g(1) = c_1,\] and such that \[d(\gamma_g(s),\gamma_g(r)) = |s-r|\,d(c_0,c_1)\quad \forall s,r\in[0,1].\] Since \(\Psi\) is a bijection, this geodesic is obtained by pulling back the straight segment between \(\Psi(c_0)\) and \(\Psi(c_1)\) in \({L^2(\Omega)}\) \[\label{eq:geodesics} \gamma_g(s) := \Psi^{-1}((1-s)\Psi(c_0) + s\,\Psi(c_1)), \quad 0\le s \le 1.\tag{17}\] In the Appendix Lemma 10 we prove that \({L^\Phi(\Omega)}\) endowed with the metric \(d\) is a complete metric space and that \(\gamma_g(s)\) given in 17 is a geodesic.

3 Auxiliary energy reformulation↩︎

3.1 The scalar case↩︎

We consider the minimization problem \[\label{eq:scalar95primal} \mathop{\mathrm{argmin}}_{u\in W^{1,\phi}_0(\Omega)} \int_{\Omega}\phi(|\nabla u|)\,\mathrm{d}x- F(u),\tag{18}\] where \(W^{1,\phi}_0(\Omega)\) is the closure of \(C^\infty_c(\Omega)\) with respect to the \(W^{1,\phi}(\Omega)\)-norm and \(F\) is a continuous linear operator in \(W^{1,\phi}_0(\Omega)^*\). The Euler–Lagrange equations associated with 18 read: Find \(u\in W^{1,\phi}_0(\Omega)\) such that \[\label{eq:scalar95EL} \int_{\Omega} \frac{\phi'(|\nabla u|)}{|\nabla u|}\nabla u \cdot \nabla v \,\mathrm{d}x = F(v) \quad \forall v\in W^{1,\phi}_0(\Omega).\tag{19}\] For regularity results of the minimizers of 18 under Assumption 4, see [49], [50] the references therein.

Using the definition of \(\Phi\) from 7 , we rewrite the energy in the equivalent form \[\label{eq:scalar95primal95Phi} E(u):=\int_{\Omega}\Phi(|\nabla u|^2)\,\mathrm{d}x- F(u),\tag{20}\] yielding the Euler–Lagrange equations: Find \(u\in W^{1,\phi}_0(\Omega)\) such that \[\label{eq:scalar95EL95mu} \int_{\Omega}\mu(|\nabla u|^2)\,\nabla u\cdot\nabla v\,\mathrm{d}x = F(v) \quad \forall v\in W^{1,\phi}_0(\Omega).\tag{21}\]

As stated in the introduction, our goal is not to derive an equivalent minimization problem, but to introduce an auxiliary-field formulation that shifts the nonlinearity from the state equation to a scalar variable, and that provides a variationally consistent framework for designing and supporting approximation schemes through a sequence of uniformly elliptic linear state problems. To this end, we first introduce the extended functional \[J:{W_0^{1,\phi}(\Omega)}\times {L^\Phi(\Omega)}\rightarrow [-\infty,+\infty], \quad J(u,c) := \begin{cases} \displaystyle \int_{\Omega} \Phi(c)+\frac{1}{2}\,\mu(c)\big(|\nabla u|^2-c\big)\,\mathrm{d}x- F(u), & c\in{L_0^\Phi(\Omega)},\\ +\infty,&\text{otherwise}, \end{cases}\] where \[\label{eq:L95Phi} {L_0^\Phi(\Omega)}:= \left\{ c \in {L^\Phi(\Omega)}: c(x) \ge 0 \,\text{ a.e. in }\,\Omega \right\}.\tag{22}\] For a fixed \(u\in {W_0^{1,\phi}(\Omega)}\), since \(2\Phi'(c) = \mu(c)\), stationarity of \(J\) with respect to \(c\) implies \[\mu'(c)(|\nabla u|^2-c) = 0,\quad\text{ a.e. in }\Omega,\] and thus, from Assumption 3, we obtain \[|\nabla u|^2 = c\, \text{ a.e. in }\{x\in\Omega :c(x) > 0\}.\] Stationarity of \(J\) with respect to \(u\) instead yields \[\label{eq:scalar95state95nwp} a[c](u,v) = F(v), \quad \forall v\in {W_0^{1,\phi}(\Omega)},\tag{23}\] where \[\label{eq:scalar95bilin} a[c](u,v):=\int_{\Omega}\mu(c)\,\nabla u\cdot\nabla v\,\mathrm{d}x.\tag{24}\] For a general \(c\in {L_0^\Phi(\Omega)}\), the variational problem 23 is not well-posed on \({W_0^{1,\phi}(\Omega)}\) and results like the Minty-Browder theorem [51] cannot be applied. In the concave case we can in fact prove the coercivity of the functional in \({W_0^{1,\phi}(\Omega)}\), but the bilinear form 24 is not well defined since \(\mu(c)\nabla\cdot\) does not map \({W_0^{1,\phi}(\Omega)}\) into \({W_0^{1,\phi}(\Omega)}^*\). On the other hand, when \(\Phi\) is convex, the linear operator induced by the bilinear form is hemicontinuous, but we lose coercivity instead.

When \(\Phi\) is concave, we have \[\Phi(r) \leq \Phi(c) + \frac{1}{2} \mu(c) (r - c)\quad \forall\, r,c \ge 0.\] Integrating on \(\Omega\) and subtracting \(F(u)\) from both sides yields \[E(u) \le J(u,c), \quad \forall u\in {W_0^{1,\phi}(\Omega)},\, \forall c \in {L_0^\Phi(\Omega)}.\] and thus \[\label{eq:Eu95bound95conc} \inf_{u\in {W_0^{1,\phi}(\Omega)}}E(u) \le \inf_{u\in {W_0^{1,\phi}(\Omega)}} J(u,c) = \mathcal{J}_c(c) + \inf_{u\in {W_0^{1,\phi}(\Omega)}} \Big( \int_\Omega \frac{1}{2}\mu(c)|\nabla u|^2 \,\mathrm{d}x- F(u)\Big),\tag{25}\] where the functional \[\label{eq:Jc95def} \mathcal{J}_c(c) := \int_\Omega \Phi(c) -\frac{1}{2} \mu(c) c \,\mathrm{d}x,\tag{26}\] is well defined and finite for all \(c\in {L_0^\Phi(\Omega)}\) by Lemma 3. In the case of convex \(\Phi\), the comparison with the tangent functional reverses \[\Phi(r) \geq \Phi(c) + \frac{1}{2} \mu(c) (r - c)\quad \forall\, r,c \ge 0,\] yielding \[E(u) \ge J(u,c), \quad \forall u\in {W_0^{1,\phi}(\Omega)},\, \forall c \in {L_0^\Phi(\Omega)},\] and \[\label{eq:Eu95bound95conv} \inf_{u\in {W_0^{1,\phi}(\Omega)}}E(u) \ge \inf_{u\in {W_0^{1,\phi}(\Omega)}} J(u,c) = \mathcal{J}_c(c) + \inf_{u\in {W_0^{1,\phi}(\Omega)}} \Big( \int_\Omega \frac{1}{2}\mu(c)|\nabla u|^2 \,\mathrm{d}x- F(u)\Big).\tag{27}\]

We then define an auxiliary energy in \(c\) by replacing the space \({W_0^{1,\phi}(\Omega)}\) in the minimization of the Dirichlet integral with \({H_0^1(\Omega)}\) and constraining \(c\) to the set \[\label{eq:C95d95def} {L_\epsilon^\Phi(\Omega)}:= \left\{c \in {L^\Phi(\Omega)}:\epsilon \le c(x) \le 1/\epsilon \,\text{ a.e. in } \Omega \right\},\tag{28}\] for an arbitrary \(\epsilon \in (0,1)\), which ensures the uniform bounds to hold almost everywhere \[\label{eq:mu95unif95bounds} 0 < \mu_- \le \mu(c(x)) \le \mu_+ < +\infty, \quad \mu_- :=\min\{\mu(\epsilon), \mu(1/\epsilon)\}, \quad \mu_+ :=\max\{\mu(\epsilon), \mu(1/\epsilon)\}.\tag{29}\] For the problem to be well-posed, we need a further assumption.

Note that the above condition is always satisfied in the concave case but not in the convex case. We also point out that, for ease of discussion, from now on, all the inequalities (and equalities) involving pointwise bounds for functions must be understood in the almost everywhere sense.

Under Assumption 6 we can consider the well-posed functional for a given \(c\in{L_\epsilon^\Phi(\Omega)}\), see e.g [52], \[\label{eq:Ju95def} \mathcal{J}_u(c) := \inf_{u\in {H_0^1(\Omega)}} \left( \int_{\Omega}\frac{1}{2}\mu(c)|\nabla u|^2\,\mathrm{d}x- F(u)\right),\tag{30}\] whose Euler–Lagrange equations lead to the weighted Poisson problem: Find \(u[c] \in {H_0^1(\Omega)}\) such that \[\label{eq:scalar95state} \int_\Omega\mu(c)\nabla u[c]\cdot\nabla v \,\mathrm{d}x= F(v), \quad \forall v\in {H_0^1(\Omega)},\tag{31}\] yielding the identities \[\label{eq:energy95u95c} \int_{\Omega}\mu(c)|\nabla u[c]|^2\,\mathrm{d}x= F(u[c]),\quad \mathcal{J}_u(c) = -\frac{1}{2}F(u[c]) = -\frac{1}{2}\int_{\Omega}\mu(c)|\nabla u[c]|^2\,\mathrm{d}x.\tag{32}\] Combining 25 with the embedding 6 in the concave case we find \[\inf_{u\in {W_0^{1,\phi}(\Omega)}}E(u) \le \mathcal{J}_c(c) + \mathcal{J}_u(c),\quad \forall c\in {L_\epsilon^\Phi(\Omega)}.\] Similarly, for the convex case, we combine 27 with the embedding 5 and obtain \[\inf_{u\in {W_0^{1,\phi}(\Omega)}}E(u) \ge \mathcal{J}_c(c) +\mathcal{J}_u(c),\quad \forall c\in {L_\epsilon^\Phi(\Omega)}.\] Using the definition of \(\omega\) given in 14 we arrive at the extended auxiliary energy functional \[\label{eq:E95c} {\mathcal{E}_\epsilon}: {L^\Phi(\Omega)}\to (-\infty,+\infty], \quad {\mathcal{E}_\epsilon}(c) := \begin{cases} -\omega \mathcal{J}_c(c) -\omega \mathcal{J}_u(c),& c\in{L_\epsilon^\Phi(\Omega)},\\ +\infty,&\text{otherwise}, \end{cases}\tag{33}\] which is a proper, well-defined functional, and it is bounded from below by the value of the energy of the minimizer of 20 multiplied by \(-\omega\), since we have \[\label{eq:Ee95bounded} -\omega E(u[c]) \le -\omega J(u[c], c) ={\mathcal{E}_\epsilon}(c),\quad c \in {L_\epsilon^\Phi(\Omega)}.\tag{34}\]

3.2 The incompressible vector-valued case↩︎

The extension of the previous results to the vector-valued case \[\mathop{\mathrm{argmin}}_{{\boldsymbol{u}}\in W^{1,\phi}_{0}(\Omega;\mathbb{R}^d)}\int_{\Omega}\phi(|\nabla{\boldsymbol{u}}|)\,\mathrm{d}x- F({\boldsymbol{u}}),\] is straightforward. In this section, we instead consider the minimization problem \[\label{eq:vector95min} \mathop{\mathrm{argmin}}_{{\boldsymbol{u}}\in W^{1,\phi}_{0,\mathrm{div}}(\Omega;\mathbb{R}^d)}\int_{\Omega}\phi(|\varepsilon({\boldsymbol{u}})|)\,\mathrm{d}x- F({\boldsymbol{u}}),\tag{35}\] where \[W^{1,\phi}_{0,\mathrm{div}}(\Omega;\mathbb{R}^d) :=\{\boldsymbol{v}\in {W_0^{1,\phi}(\Omega;\mathbb{R}^d)}: \nabla\cdot\boldsymbol{v}=0\}, \quad \varepsilon({\boldsymbol{u}}):=\frac{1}{2}(\nabla{\boldsymbol{u}}+\nabla{\boldsymbol{u}}^{\mathsf T}),\] and \(F\in W^{1,\phi}_0(\Omega;\mathbb{R}^d)^\ast\). For existence, uniqueness, and regularity of the solutions under Assumption 4 see the recent work [53] and the references therein.

The auxiliary functional is constructed by following verbatim the scalar case, replacing \(|\nabla u|^2\) by \(|\varepsilon({\boldsymbol{u}})|^2\), and the infimum of the variational problem \[\inf_{{\boldsymbol{u}}\in W^{1,\phi}_{0,\mathrm{div}}(\Omega;\mathbb{R}^d)} \left(\int_\Omega \frac{1}{2}\mu(c)|\varepsilon({\boldsymbol{u}})|^2 \,\mathrm{d}x-F({\boldsymbol{u}})\right),\] is replaced by \[\label{eq:Ju95def95vector} \mathcal{J}_{\boldsymbol{u}}(c) := \inf_{{\boldsymbol{u}}\in H^{1}_{0,\mathrm{div}}(\Omega;\mathbb{R}^d)} \left(\int_\Omega \frac{1}{2}\mu(c)|\varepsilon({\boldsymbol{u}})|^2 \,\mathrm{d}x-F({\boldsymbol{u}})\right),\tag{36}\] for a given \(c\in{L_\epsilon^\Phi(\Omega)}\) under the following assumption:

To model incompressibility, we introduce a Lagrange multiplier through the pressure variable that lives in \[Q:=\left\{p\in {L^2(\Omega)}: \int_\Omega p \,\mathrm{d}x= 0\right\},\] and obtain the Euler–Lagrange equations for 36 given by the weighted incompressible Stokes system: Find \({\boldsymbol{u}}[c] \in {H_0^1(\Omega;\mathbb{R}^d)}\) and \(p[c] \in Q\) such that \[\begin{align} \boldsymbol{a}[c]({\boldsymbol{u}}[c], \boldsymbol{v})- b(\boldsymbol{v}, p[c]) &= F(\boldsymbol{v}), \quad \forall \boldsymbol{v}\in {H_0^1(\Omega;\mathbb{R}^d)}, \tag{37} \\ b({\boldsymbol{u}}[c],q) &= 0, \quad \forall q\in Q, \tag{38} \end{align}\] where \[\label{eq:vector95bilin} \boldsymbol{a}[c]({\boldsymbol{u}}, \boldsymbol{v}) := \int_{\Omega}\mu(c)\,\varepsilon({\boldsymbol{u}}):\varepsilon(\boldsymbol{v})\,\mathrm{d}x, \quad b({\boldsymbol{u}},q) := \int_{\Omega} q\,\nabla\cdot{\boldsymbol{u}}\,\mathrm{d}x.\tag{39}\] The well-posedness of the weighted incompressible Stokes problem is standard for a given \(c\in{L_\epsilon^\Phi(\Omega)}\), see e.g [52]. Testing 37 with \(\boldsymbol{v}={\boldsymbol{u}}[c]\) and using the incompressibility of \({\boldsymbol{u}}[c]\) gives the identities \[\label{eq:energy95u95c95vector} \int_{\Omega}\mu(c)|\varepsilon({\boldsymbol{u}}[c])|^2\,\mathrm{d}x= F({\boldsymbol{u}}[c]),\quad \mathcal{J}_{\boldsymbol{u}}(c) = -\frac{1}{2}F({\boldsymbol{u}}[c]) = -\frac{1}{2}\int_{\Omega}\mu(c)|\varepsilon({\boldsymbol{u}}[c])|^2\,\mathrm{d}x,\tag{40}\] and the auxiliary functional has the same structure as the scalar case \[\label{eq:E95c95vector} {\mathcal{E}_\epsilon}: {L^\Phi(\Omega)}\to (-\infty,+\infty], \quad {\mathcal{E}_\epsilon}(c) := \begin{cases} -\omega \mathcal{J}_c(c) -\omega \mathcal{J}_{\boldsymbol{u}}(c),& c\in{L_\epsilon^\Phi(\Omega)},\\ +\infty,&\text{otherwise}. \end{cases}\tag{41}\] Even in this case, \({\mathcal{E}_\epsilon}\) is bounded from below by the minimum of 35 multiplied by \(-\omega\) and we have \(-\omega E({\boldsymbol{u}}[c]) \le {\mathcal{E}_\epsilon}(c)\) for all \(c \in {L_\epsilon^\Phi(\Omega)}\). Stationarity with respect to a given \(c\) yields \[\mu'(c)(|\varepsilon({\boldsymbol{u}})|^2-c)=0 \quad \text{a.e. in }\Omega.\]

3.3 Lower semicontinuity of the auxiliary functional↩︎

We then prove the lower semicontinuity of \({\mathcal{E}_\epsilon}\). First we prove the continuity of \(\mathcal{J}_c\) on \({L_0^\Phi(\Omega)}\); then, we prove the continuity of \(\mathcal{J}_u\) and \(\mathcal{J}_{\boldsymbol{u}}\) on \({L_\epsilon^\Phi(\Omega)}\). Note that both spaces \({L_0^\Phi(\Omega)}\) and \({L_\epsilon^\Phi(\Omega)}\) are closed in the \(d-\)topology.

We use the notation \(h(r) := \Phi(r) -\frac{1}{2}\mu(r)r\), \(h'(r) = - \frac{1}{2}\mu'(r) r\); hence \[h(r) = \frac{1}{2}\int^r_0 -\mu'(s)s\,\mathrm{d}s.\] Combining 10 with Assumption 4 yields \[\label{eq:temp95muprime95bound} A \mu(s) \le \omega s \mu'(s) \le B \mu(s),\tag{42}\] with \(A = \frac{{{S}^{-}_{\phi}}-1}{2}\) and \(B = \frac{{{S}^{+}_{\phi}}-1}{2}\) when \(\Phi\) is convex, while \(A = \frac{1-{{S}^{+}_{\phi}}}{2}\) and \(B = \frac{1-{{S}^{-}_{\phi}}}{2}\) when \(\Phi\) is concave. Integrating both sides from zero to \(r\), we obtain \(A \Phi(r) \le -\omega h(r) \le B \Phi(r)\), and thus \(-\omega \mathcal{J}_c(c)\) is positive and controlled by the modular of \(c\) in \({L_0^\Phi(\Omega)}\).

Defining \(q(\eta) := -\omega h(\psi^{-1}(\eta))\) with \(\eta := \psi(r)\) and \(\psi\) given by 13 , straightforward calculations yield \[\label{eq:k95eta95prime95temp} \begin{align} q'(\eta) = -\omega \frac{h'(r)}{\psi'(r)} = r\frac{\omega \mu'(r)/2}{\sqrt{\omega \mu'(r)/2}} = \sqrt{\frac{r}{2}}\sqrt{\omega \mu'(r)r} \leq \sqrt{\frac{B}{2}}\sqrt{\mu(r) r} \leq L \psi(r) = L\eta, \end{align}\tag{43}\] where the first inequality follows from the upper bound in 42 while the second inequality follows from the lower bound in statement (i) of Lemma [lem:psi95bounds95mu], thus the constant is \(L=\sqrt{\frac{B}{2C_1}}\). From 43 , we obtain \[\label{eq:k95eta95prime95temp952} \left|q(\eta_2)-q(\eta_1)\right| = \left|\int^{\eta_2}_{\eta_1}q'(s)\,\mathrm{d}s\right|\leq \frac{L}{2}\left|\eta^2_2 - \eta^2_1\right| = \frac{L}{2}\left|\eta_2 - \eta_1\right|\left|\eta_2 + \eta_1\right|,\quad \forall \eta_1,\eta_2 \ge 0.\tag{44}\] Let \(\{c_n\}_{n\in\mathbb{N}}\in {L_0^\Phi(\Omega)}\) such that \(c_n\to c^*\) in the \(d\)-topology; then, since \({L_0^\Phi(\Omega)}\) is closed, \(c^\ast \in {L_0^\Phi(\Omega)}\). Using the notation \(\eta_n := \Psi(c_n)\), \(\eta^* := \Psi(c^*)\), we have, using 44 , the Cauchy–Schwarz inequality, a triangle inequality, and the definition of the metric 16 \[\begin{align} |\mathcal{J}_c(c_n) - \mathcal{J}_c(c^*)| &\leq \frac{L}{2} \int_\Omega \left|\eta_n - \eta^*\right| \left|\eta_n + \eta^*\right| \,\mathrm{d}x\\ &\leq \frac{L}{2} \|\eta_n +\eta^*\|_{L^2(\Omega)}\|\eta_n -\eta^*\|_{L^2(\Omega)}\leq \frac{L}{2} (d(c_n,c^*) + 2\|\eta^*\|_{L^2(\Omega)}) d(c_n,c^*) \to 0, \end{align}\] and thus \(-\omega \mathcal{J}_c\) is continuous on \({L_0^\Phi(\Omega)}\).

Let \(\{c_n\}_{n\in\mathbb{N}}\in {L_\epsilon^\Phi(\Omega)}\) such that \(c_n\to c^*\) in the \(d\)-topology; then, since \({L_\epsilon^\Phi(\Omega)}\) is closed in the \(d\)-topology, \(c^\ast \in {L_\epsilon^\Phi(\Omega)}\). By Chebyshev’s inequality, \(\Psi(c_n) \to \Psi(c^\ast)\) in measure on \(\Omega\). Since \(\Psi^{-1}\) is continuous, we also have \(c_n \to c^\ast\) in measure, as well as \(\mu(c_n) \to \mu(c^\ast)\) since \(\mu\) continuous and finite on \({L_\epsilon^\Phi(\Omega)}\) for 29 .

We first give the proof for the scalar case, and show that \(u[c_n]\) converges to \(u[c^\ast]\) in \({H_0^1(\Omega)}\); this allows us to prove the statement since we obtain \(F(u[c_n])\to F(u[c^\ast])\) given that \(F\in{H_0^1(\Omega)}^\ast\) and, using the identity in 32 , \(-\omega \mathcal{J}_u(c_n)\to -\omega \mathcal{J}_u(c^\ast)\). For convenience, we write \(a_n:=\mu(c_n)\), \(a^\ast:=\mu(c^\ast)\), \(u_n:=u[c_n]\), \(u^\ast:=u[c^\ast]\), then, from 31 \[\int_\Omega a_n\nabla u_n\cdot\nabla v\,\mathrm{d}x=F(v)=\int_\Omega a^\ast\nabla u^\ast\cdot\nabla v \,\mathrm{d}x, \quad\forall v\in{H_0^1(\Omega)}.\] Subtracting \(\int_\Omega a_n\nabla u^\ast\cdot\nabla v\,\mathrm{d}x\) from both sides yields \[\int_\Omega a_n\nabla(u_n-u^\ast)\cdot\nabla v\,\mathrm{d}x=\int_\Omega (a^\ast-a_n)\nabla u^\ast\cdot\nabla v\,\mathrm{d}x.\] Choosing \(v=u_n-u^\ast\) we obtain \[\int_\Omega a_n|{\nabla(u_n-u^\ast)}|^2\,\mathrm{d}x=\int_\Omega (a^\ast-a_n)\nabla u^\ast\cdot\nabla(u_n-u^\ast)\,\mathrm{d}x.\] Using the lower bound \(a_n\ge \mu_-\) from 29 , and a Cauchy–Schwarz inequality argument we get \[\|\nabla(u_n-u^\ast)\|_{L^2(\Omega)}\le\frac{1}{\mu_-}\|(a_n-a^\ast)\nabla u^\ast\|_{L^2(\Omega)}.\] The right-hand side of the above inequality is uniformly integrable since \(u^\ast\in {H_0^1(\Omega)}\) and \[|(a_n-a^\ast)\nabla u^\ast|^2\le 4\mu_+^2|\nabla u^\ast|^2.\] Therefore, since \(a_n\to a^\ast\) in measure, by Lebesgue-Vitali theorem [54], \(\|\nabla(u_n-u^\ast)\|_{L^2(\Omega)}\to 0\). By Poincaré’s inequality, convergence of the gradients implies \(u[c_n]\to u[c^\ast]\) in \({H_0^1(\Omega)}\).

The proof for the incompressible vector-valued case is very similar. Using the shorthand notation \({\boldsymbol{u}}_n:={\boldsymbol{u}}[c_n]\), \({\boldsymbol{u}}^\ast:={\boldsymbol{u}}[c^\ast]\), \(p_n:=p[c_n]\), and \(p^\ast:=p[c^\ast]\), from 37 and 38 we obtain \[\begin{align} \int_\Omega a_n\varepsilon({\boldsymbol{u}}_n-{\boldsymbol{u}}^\ast)\!:\!\varepsilon(\boldsymbol{v})\,\mathrm{d}x-\int_\Omega\nabla\cdot\boldsymbol{v}\,(p_n - p^\ast)\,\mathrm{d}x&=\int_\Omega (a^\ast-a_n)\,\varepsilon({\boldsymbol{u}}^\ast)\!:\!\varepsilon(\boldsymbol{v})\,\mathrm{d}x,\quad \forall\boldsymbol{v}\in{H_0^1(\Omega;\mathbb{R}^d)},\\ \int_\Omega\nabla\cdot({\boldsymbol{u}}_n - {\boldsymbol{u}}^\ast)\,q\,\mathrm{d}x&= 0, \quad \forall q\in Q. \end{align}\] Choosing the divergence-free test function \(\boldsymbol{v}= {\boldsymbol{u}}_n - {\boldsymbol{u}}^\ast\) yields \[\int_\Omega a_n|{\varepsilon({\boldsymbol{u}}_n-{\boldsymbol{u}}^\ast)}|^2\,\mathrm{d}x=\int_\Omega (a^\ast-a_n)\varepsilon({\boldsymbol{u}}^\ast)\!:\!\varepsilon({\boldsymbol{u}}_n-{\boldsymbol{u}}^\ast)\,\mathrm{d}x.\] Repeating the same arguments, using the symmetric gradient instead of the plain gradient, we arrive at \(\|\varepsilon({\boldsymbol{u}}_n-{\boldsymbol{u}}^\ast)\|_{L^2(\Omega)}\to 0\) and thus, by Korn’s and Poincaré’s inequalities, \({\boldsymbol{u}}[c_n]\to {\boldsymbol{u}}[c^\ast]\) in \({H_0^1(\Omega;\mathbb{R}^d)}\).

Lower semicontinuity accounts for proving that \[{\mathcal{E}_\epsilon}(c^\ast)\le \liminf_{n\to+\infty}{\mathcal{E}_\epsilon}(c_n)\] for all sequences \(\{c_n\}_{n\in\mathbb{N}}\in {L^\Phi(\Omega)}\) that converge to \(c^\ast \in {L^\Phi(\Omega)}\). We split the proof into two cases.

Suppose \(c^\ast\notin {L_\epsilon^\Phi(\Omega)}\), then \({\mathcal{E}_\epsilon}(c^\ast)=+\infty\). Since \({L_\epsilon^\Phi(\Omega)}\) is closed with respect to \(d\), then eventually \(c_n\notin {L_\epsilon^\Phi(\Omega)}\) because the complement of \({L_\epsilon^\Phi(\Omega)}\) is open, and thus \(\liminf_{n\to+\infty}{\mathcal{E}_\epsilon}(c_n)=+\infty\).

Suppose now that \(c^\ast\in {L_\epsilon^\Phi(\Omega)}\). If \(c_n\in{L_\epsilon^\Phi(\Omega)}\) only for finitely many terms, then \(\liminf_{n\to+\infty}{\mathcal{E}_\epsilon}(c_n)=+\infty\), and thus, for such sequences, the functional is lower semicontinuous. It only remains to consider the case where \(c_n\in{L_\epsilon^\Phi(\Omega)}\) for infinitely many terms, having indices \(I:=\left\{n_k \in \mathbb{N} :c_{n_k} \in {L_\epsilon^\Phi(\Omega)}\right\}\). Since the terms of the sequences for which \(n\notin I\) do not lower the limit inferior, we have \[\liminf_{n\to+\infty}{\mathcal{E}_\epsilon}(c_n) = \liminf_{k\to+\infty}{\mathcal{E}_\epsilon}(c_{n_k}).\] From Theorems 4 and 5, \({\mathcal{E}_\epsilon}\) is continuous on \({L_\epsilon^\Phi(\Omega)}\), and thus \({\mathcal{E}_\epsilon}(c_{n_k})\to {\mathcal{E}_\epsilon}(c^\ast)\), proving that the functional is lower semicontinuous.

3.4 Gradient flow in metric space↩︎

In this section, we study the minimization of the auxiliary energy \({\mathcal{E}_\epsilon}\) using the theory of gradient flows in metric spaces [35]. Using the notation \(C:={L^\Phi(\Omega)}\), we consider the metric functional system as in [55] \[\label{eq:metsym} \begin{align} &(C,d,{\mathcal{E}_\epsilon}), \text{ where } (C,d) \text{ is a complete metric space},\\ &{\mathcal{E}_\epsilon}: C\rightarrow (-\infty,+\infty]\text{ is a proper, lower semicontinuous functional, bounded from below}, \end{align}\tag{45}\] and provide sufficient conditions for the existence of Evolution Variational Inequalities (EVI). We first recall the definition of \(\mathrm{EVI}_\lambda\)-gradient flows.

We recall here some of the properties of the solutions \(\mathrm{EVI}_\lambda\)-gradient flows taken from [56] when the space \((C,d)\) is complete \[\label{eq:evil95props} \begin{align} \boldsymbol{contraction property}:\quad& d(S_t(c_0),S_t(c_1))\leq e^{-\lambda(t-s)}d(S_s(c_0),S_s(c_1)),\quad 0\leq s\leq t < +\infty,\\ \boldsymbol{asymptotic behaviour}:\quad& \text{If }\lambda > 0 \Rightarrow d(S_t(c),c^*) \leq e^{-\lambda(t-t_0)}d(S_{t_0}(c),c^*),\quad c^\ast = \mathop{\mathrm{argmin}}_{c\in C}{\mathcal{E}_\epsilon}(c). \end{align}\tag{46}\] The contraction estimate gives uniqueness and continuous dependence on the initial datum. If \(\lambda\ge 0\), the distance between two solutions is nonincreasing, while if \(\lambda>0\), it decays exponentially. In the latter case, the asymptotic estimate shows that the flow selects the unique minimizer \(c^*\) of \({\mathcal{E}_\epsilon}\).

The evolution maps associated with the \(\mathrm{EVI}_\lambda\)-gradient flows are approximated by interpolants obtained by incremental minimizations using the method of minimizing movements [36]; at iteration \(k\), we define iteratively \[\label{eq:mm} c_{k+1} \in \min_{c\in C} \,{\mathcal{E}_\epsilon}(c) + \frac{1}{2\tau_k} d(c_k, c)^2,\tag{47}\] where \(\tau_k\) is the current time step, and we define the left-continuous piecewise constant interpolant \[\overline{c}_{\boldsymbol{\tau}}(t) := \begin{cases} c_0, &\quad t=0,\\ c_{k+1} &\quad t \in(t_k, t_{k+1}], \end{cases}\] where discrete times are given recursively by \(t_{k+1} = t_k + \tau_k\), with \(t_0 = 0\). Under the existence of an \(\mathrm{EVI}_\lambda\) flow with \(\lambda > 0\), the following a priori error estimate holds [35] for a fixed time step \(\tau\) \[d(\overline{c}_{\boldsymbol{\tau}}(t), S_t(c_0)) \leq K(c_0) \sqrt{\tau} e^{-\lambda_\tau t}, \quad \lambda_\tau := \frac{\log (1+\lambda \tau)}{\tau},\] where \(K(c_0)\) is a constant depending on \(c_0\). Combining the estimate above with the asymptotic behavior of \(\mathrm{EVI}_\lambda\) flows from 46 gives \[\label{eq:apriori95lambda95asym} d(\overline{c}_{\boldsymbol{\tau}}(t), c^\ast) \leq d(\overline{c}_{\boldsymbol{\tau}}(t), S_t(c_0)) + d(S_t(c_0), c^*)\leq K(c_0) \sqrt{\tau} e^{-\lambda_\tau t} + d(c_0, c^*)e^{-\lambda t}\le Ke^{-\lambda_\tau t}.\tag{48}\]

Central to the proof of the existence of the solutions of evolutionary variational inequalities is the notion of geodesical \(\lambda\)-convexity of a functional, which is stated below.

Under the conditions that the space \((C,d)\) is a complete geodesic space and that the map \(c\mapsto\frac{1}{2}d(\bar{c},c)^2\) is geodesically \(1\)-convex for all \(\bar{c}\in C\), proving the existence of \(\mathrm{EVI}_\lambda\)-gradient flows reduce to prove the geodesical \(\lambda\)-convexity of \({\mathcal{E}_\epsilon}\) [55]–[57]. The geodesical 1-convexity of the metric 16 follows from an application of the Hilbert-space identity in \({L^2(\Omega)}\) for the squared distance transported through \(\Psi\). The following Theorem proves the geodesical \(\lambda\)-convexity of \({\mathcal{E}_\epsilon}\) under sufficient assumptions.

Given that \({\mathcal{E}_\epsilon}=+\infty\) on \({L_\epsilon^\Phi(\Omega)}^c\), we only need to prove geodesical \(\lambda\)-convexity of \({\mathcal{E}_\epsilon}\) on \({L_\epsilon^\Phi(\Omega)}\); instead of proving it directly on \({L_\epsilon^\Phi(\Omega)}\), we leverage the isometry of \({L^\Phi(\Omega)}\) and \({L^2(\Omega)}\) via the map \(\Psi\) given in 15 and prove the equivalent statement that the functional \[\label{eq:E95l2} \widetilde{{\mathcal{E}_\epsilon}}: {L^2(\Omega)}\longrightarrow(-\infty,+\infty], \quad \widetilde{{\mathcal{E}_\epsilon}}(\eta):=({\mathcal{E}_\epsilon}\circ\Psi^{-1})(\eta) ,\tag{49}\] is \(\lambda\)-convex on the image of \({L_\epsilon^\Phi(\Omega)}\) under \(\Psi\), that is the closed convex interval \[{L_\epsilon^2(\Omega)}:=\left\{\eta \in {L^2(\Omega)}:\psi(\epsilon)\le \eta\le \psi(1/\epsilon)\right\},\] where \(\psi\) is given in 13 . Specifically, we derive sufficient conditions under which the function \[H(t):=\widetilde{{\mathcal{E}_\epsilon}}(\eta_t)-\frac{\lambda}{2}\|\eta_t\|_{L^2(\Omega)}^2,\quad \eta_t := (1-t) \eta_0 + t\eta_1,\] is convex, which allows us to prove the \(\lambda\)-convexity of \(\widetilde{{\mathcal{E}_\epsilon}}\) in \({L_\epsilon^2(\Omega)}\). In fact, if \[H(t) \le (1-t) H(0) + tH(1) = (1-t) \widetilde{{\mathcal{E}_\epsilon}}(\eta_0) + t\widetilde{{\mathcal{E}_\epsilon}}(\eta_1) -\frac{\lambda}{2}\left((1-t)\|\eta_0\|_{L^2(\Omega)}^2 + t\|\eta_1\|_{L^2(\Omega)}^2\right),\] using the identity \[\|\eta_t\|_{L^2(\Omega)}^2 = (1-t)\|\eta_0\|_{L^2(\Omega)}^2 + t\|\eta_1\|_{L^2(\Omega)}^2 -t(1-t)\|\eta_1-\eta_0\|_{L^2(\Omega)}^2,\] yields \[\widetilde{{\mathcal{E}_\epsilon}}(\eta_t) \le (1-t) \widetilde{{\mathcal{E}_\epsilon}}(\eta_0) + t\widetilde{{\mathcal{E}_\epsilon}}(\eta_1)-\frac{\lambda}{2}t(1-t)\|\eta_1 - \eta_0\|_{L^2(\Omega)}^2,\] and thus the geodesical \(\lambda\)-convexity of \({\mathcal{E}_\epsilon}\) \[{\mathcal{E}_\epsilon}(c_t) \le (1-t) {\mathcal{E}_\epsilon}(c_0) + t{\mathcal{E}_\epsilon}(c_1)-\frac{\lambda}{2}t(1-t)d(c_1,c_0)^2,\] where \(c_t := \Psi^{-1}(\eta_t)\) is the geodesic in \({L_\epsilon^\Phi(\Omega)}\).

In the sequel, we use the notation \(K_\epsilon:=[\epsilon,1/\epsilon]\), \(I_\epsilon:=[\psi(\epsilon), \psi(1/\epsilon)]\), together with \(\xi:=\eta_1-\eta_0\). Note that \(\xi\in{L^\infty(\Omega)}\); moreover, since \(\mu\in C^2(0,+\infty)\) and \(\omega\mu'\) is strictly positive on \(K_\epsilon\), we have that \(\psi\) is a \(C^2\) diffeomorphisms from \(K_\epsilon\) onto \(I_\epsilon\).

We first consider the functional \[\widetilde{\mathcal{J}}_c:=-\omega(\mathcal{J}_c\circ\Psi^{-1})(\eta),\] and define, as in the proof of Theorem 4, \[h(r) := \Phi(r) -\frac{1}{2}\mu(r)r,\quad q(s):=-\omega h(\psi^{-1}(s)), \quad s := \psi(r)\,\] so that \[\widetilde{\mathcal{J}}_c(\eta)=\int_\Omega q(\eta)\,\mathrm{d}x.\] Since \(\eta_t\) takes values in \(I_\epsilon\), and \(q\in C^2(I_\epsilon)\), the chain rule gives \[\label{eq:J95c952} \frac{d^2}{dt^2} \widetilde{\mathcal{J}}_c(\eta_t)=\int_\Omega q''(\eta_t)\xi^2\,\mathrm{d}x.\tag{50}\] It remains to compute \(q''\): since \(\displaystyle \frac{dr}{ds} = \frac{1}{\sqrt{\omega\mu'(r)/2}}\), we have \[\begin{align} q'(s) &= -\omega h'(r) \frac{dr}{ds} =\frac{1}{2}\omega \mu'(r)r \frac{1}{\sqrt{\omega\mu'(r)/2}} = r\sqrt{\omega\mu'(r)/2},\nonumber\\ q''(s) &= \frac{dr}{ds}\sqrt{\omega\mu'(r)/2} + r \frac{1}{2\sqrt{\omega\mu'(r)/2}}\frac{\omega \mu''(r)}{2}\frac{dr}{ds} = 1 +\frac{r\mu''(r)}{2\mu'(r)}\label{eq:temp95q952}. \end{align}\tag{51}\]

We then consider the term \[\label{eq:Jueta} \widetilde{\mathcal{J}}_u(\eta) = \frac{\omega}{2}\int_\Omega \widetilde{\mu}(\eta)|\nabla u[\eta]|^2\,\mathrm{d}x, \quad \widetilde{\mu}(\eta) := (\mu \circ\Psi^{-1})(\eta),\tag{52}\] under the constraint \[\label{eq:Jueta95state} \int_\Omega \widetilde{\mu}(\eta)\nabla u[\eta]\cdot\nabla v\,\mathrm{d}x= F(v),\quad \forall v \in {H_0^1(\Omega)}.\tag{53}\] We first treat the scalar case in detail and outline the changes needed to extend the proof to the incompressible case later. For \(t\in[0,1]\), we first prove that the map \(t\mapsto u[\eta_t]\) is continuous from \([0,1]\) into \(V:={H_0^1(\Omega)}\) and twice differentiable from \((0,1)\) into \(V\). Since \(\xi\in L^\infty(\Omega)\) and \(\widetilde{\mu}\in C^2(I_\epsilon)\), the map \(t\mapsto \widetilde{\mu}(\eta_t)\) is \(C^2((0,1);L^\infty(\Omega))\), with \[\frac{d}{dt}\widetilde{\mu}(\eta_t)=\widetilde{\mu}'(\eta_t)\xi, \quad \frac{d^2}{dt^2}\widetilde{\mu}(\eta_t)=\widetilde{\mu}''(\eta_t)\xi^2.\] For \(t\in[0,1]\), define \[A_t:V\to V^\ast,\quad\langle A_t w,v\rangle:=\int_\Omega\widetilde{\mu}(\eta_t)\nabla w\cdot\nabla v\,\mathrm{d}x,\] with coefficient bounds \(\mu_-\le \widetilde{\mu}(\eta_t)\le\mu_+\) from 29 . Therefore \(A_t\) is an isomorphism from \(V\) onto \(V^\ast\), uniformly in \(t\). Consequently, \(t\mapsto A_t\) is \(C^2((0,1);\mathcal{L}(V,V^\ast))\). Since \(A_t\) is uniformly coercive and invertible for every \(t\), the map \(t\mapsto A_t^{-1}\) is \(C^2((0,1);L(V^\ast,V))\), and hence \(t\mapsto u[\eta_t]\in C^2((0,1);V)\).

Differentiating 52 along the geodesic gives \[\label{eq:Jueta95fvar95temp} \frac{d}{dt}\widetilde{\mathcal{J}}_u(\eta_t)=\omega\int_\Omega\widetilde{\mu}(\eta_t)\nabla u[\eta_t]\cdot\nabla\zeta_t\,\mathrm{d}x+\frac{\omega}{2}\int_\Omega\widetilde{\mu}'(\eta_t)|\nabla u[\eta_t]|^2\,\xi\,\mathrm{d}x,\tag{54}\] where we have used the notation \(\zeta_t:=\frac{d}{dt}u[\eta_t]\). Differentiating 53 yields \[\label{eq:Jueta95state95var} \int_\Omega \widetilde{\mu}(\eta_t)\nabla\zeta_t \cdot\nabla v\,\mathrm{d}x= -\int_\Omega \widetilde{\mu}'(\eta_t) \nabla u[\eta_t]\cdot \nabla v \,\xi \,\mathrm{d}x, \quad \forall v\in V.\tag{55}\] Using the above identity with \(v=u[\eta_t]\), and substituting the first term in 54 yields \[\frac{d}{dt}\widetilde{\mathcal{J}}_u(\eta_t) = -\frac{\omega}{2}\int_\Omega\widetilde{\mu}'(\eta_t)|\nabla u[\eta_t]|^2 \xi\,\mathrm{d}x.\] Differentiating once more \[\frac{d^2}{dt^2}\widetilde{\mathcal{J}}_u(\eta_t) = -\omega\int_\Omega \widetilde{\mu}'(\eta_t)\nabla u[\eta_t]\cdot \nabla\zeta_t \,\xi\,\mathrm{d}x-\frac{\omega}{2}\int_\Omega\widetilde{\mu}''(\eta_t)|\nabla u[\eta_t]|^2 \,\xi^2\,\mathrm{d}x.\] Using \(v=\zeta_t\) in 55 and substituting yields \[\label{eq:J95u952} \frac{d^2}{dt^2}\widetilde{\mathcal{J}}_u(\eta_t) = \omega\int_\Omega \widetilde{\mu}(\eta_t)|\nabla\zeta_t|^2 \,\mathrm{d}x-\frac{\omega}{2}\int_\Omega\widetilde{\mu}''(\eta_t)|\nabla u[\eta_t]|^2 \,\xi^2\,\mathrm{d}x.\tag{56}\] Using the notation \(\widetilde{\mu}(s) = \mu(r)\), with \(s=\psi(r)\), straightforward calculations show \[\begin{align} \widetilde{\mu}'(s) &= \frac{\mu'(r)}{\sqrt{\omega \mu'(r)/2}} = \omega\sqrt{2 \omega \mu'(r)},\tag{57}\\ \widetilde{\mu}''(s) &= \frac{\omega}{2\sqrt{2 \omega \mu'(r)}}2\omega \mu''(r)\frac{1}{\sqrt{\omega \mu'(r)/2}} = \omega\frac{\mu''(r)}{\mu'(r)}.\tag{58} \end{align}\]

We can now prove sufficient conditions for the convexity of \(H(t)\); from 50 and 56 we obtain \[H''(t)=\int_\Omega (q''(\eta_t) -\lambda)\xi^2\,\mathrm{d}x+ \omega\int_\Omega\widetilde{\mu}(\eta_t)|\nabla\zeta_t|^2\,\mathrm{d}x -\frac{\omega}{2}\int_\Omega\widetilde{\mu}''(\eta_t)|\nabla u[\eta_t]|^2\xi^2\,\mathrm{d}x.\] We separate the cases of convex and concave \(\Phi\). In the convex case, \(\omega=1\), and thus we can drop the second integral in the formula of \(H''(t)\) and find the pointwise sufficient condition \[\lambda \le q''(\eta_t) -\frac{1}{2}\left(\widetilde{\mu}''(\eta_t)\right)_+|\nabla u[\eta_t]|^2.\] where \(\left(x\right)_+ = \max(x,0)\). From 51 and 58 , and using the fact that \(|\nabla u[\eta_t]| =|\nabla u[c_t]|\) with \(c_t=\Psi^{-1}(\eta_t)\) the geodesic in \({L_\epsilon^\Phi(\Omega)}\), we find the sufficient condition \[\label{eq:lambda95geo95temp95conv} \lambda_{c_t} = \inf_{t\in[0,1]}\mathop{\mathrm{ess\,inf}}_{x\in\Omega}\left\{1+\frac{c_t(x)\mu''(c_t(x))}{2\mu'(c_t(x))}-\left(\frac{\mu''(c_t(x))}{2\mu'(c_t(x))}\right)_+|\nabla u[c_t](x)|^2\right\} > -\infty.\tag{59}\] Since \(\mu'(r)\ge0\) when \(\Phi\) is convex, then, if \(\mu''(r)\le 0\), we can obtain a uniform expression for \(\lambda\) that holds true for all geodesics \(c_t\) \[\lambda = \inf_{r\in [\varepsilon,1/\varepsilon]}\left\{1+\frac{r\mu''(r)}{2\mu'(r)}\right\},\] which proves statement (i) after taking the infimum for all \(\epsilon>0\). In the concave case, \(\omega=-1\), and thus we must estimate the second term in \(H''(t)\). Testing 55 with \(\zeta_t\) and using a Cauchy–Schwarz argument we obtain \[\label{eq:cs95lambdaconv95conc} \begin{align} \int_\Omega \widetilde{\mu}(\eta_t)|\nabla\zeta_t|^2\,\mathrm{d}x&= \left|\int_\Omega \widetilde{\mu}'(\eta_t) \nabla u[\eta_t]\cdot \nabla\zeta_t \,\xi \,\mathrm{d}x\right|\\ &=\left|\int_\Omega \frac{\widetilde{\mu}'(\eta_t)}{\sqrt{\widetilde{\mu}(\eta_t)}} \nabla u[\eta_t] \xi \cdot \sqrt{\widetilde{\mu}(\eta_t)}\nabla\zeta_t \,\mathrm{d}x\right|\\ &\leq\left(\int_\Omega \frac{\widetilde{\mu}'(\eta_t)^2}{\widetilde{\mu}(\eta_t)} |\nabla u[\eta_t]|^2 \xi^2\,\mathrm{d}x\right)^{1/2} \left(\int_\Omega\widetilde{\mu}(\eta_t)|\nabla\zeta_t|^2 \,\mathrm{d}x\right)^{1/2}, \end{align}\tag{60}\] and thus \[\int_\Omega \widetilde{\mu}(\eta_t)|\nabla\zeta_t|^2\,\mathrm{d}x\le \int_\Omega \frac{\widetilde{\mu}'(\eta_t)^2}{\widetilde{\mu}(\eta_t)} |\nabla u[\eta_t]|^2 \xi^2\,\mathrm{d}x.\] Plugging the above inequality into the formula for \(H''(t)\), we arrive at \[H''(t)\ge \int_\Omega (q''(\eta_t) -\lambda)\xi^2\,\mathrm{d}x+\int_\Omega \left(\frac{1}{2}\widetilde{\mu}''(\eta_t)-\frac{\widetilde{\mu}'(\eta_t)^2}{\widetilde{\mu}(\eta_t)}\right) |\nabla u[\eta_t]|^2 \xi^2\,\mathrm{d}x.\] Therefore, a sufficient pointwise condition is \[\lambda \le q''(\eta_t) - \left(\frac{\widetilde{\mu}'(\eta_t)^2}{\widetilde{\mu}(\eta_t)}-\frac{1}{2}\widetilde{\mu}''(\eta_t)\right)_+ |\nabla u[\eta_t]|^2.\] Substituting the expressions for \(q''(\eta_t)\), \(\widetilde{\mu}'(\eta_t)\) and \(\widetilde{\mu}''(\eta_t)\) from 51 , 57 and 58 we find the sufficient condition \[\label{eq:lambda95geo95temp95conc} \lambda_{c_t} = \inf_{t\in[0,1]}\mathop{\mathrm{ess\,inf}}_{x\in\Omega}\left\{1+\frac{c_t(x)\mu''(c_t(x))}{2\mu'(c_t(x))}-\left(\frac{\mu''(c_t(x))}{2\mu'(c_t(x))}-\frac{2\mu'(c_t(x))}{\mu(c_t(x))}\right)_+|\nabla u[c_t](x)|^2\right\} > -\infty.\tag{61}\] Since in the concave case \(\mu'(r)\le0\), if \(\mu''(r)\mu(r) \ge 4\mu'(r)^2\) we can obtain the simpler expression \[\lambda = \inf_{r\in [\varepsilon,1/\varepsilon]}\left\{1+\frac{r\mu''(r)}{2\mu'(r)}\right\},\] and statement (ii) is proven.

The incompressible vector-valued case is identical, up to replacing the gradient by the symmetric gradient. The differentiability of \(t\mapsto u[\eta_t]\) in the Stokes case follows equivalently, by working on the divergence-free space \(H^1_{0,\mathrm{div}}(\Omega;\mathbb{R}^d)\), where the weighted Stokes operator is uniformly coercive by Korn’s inequality and the bounds \(0<\mu_-\le \widetilde{\mu}(\eta_t)\le \mu_+\) hold. Let \(({\boldsymbol{u}}[\eta_t],p[\eta_t])\in H^1_0(\Omega;\mathbb{R}^d)\times Q\) solve the weighted Stokes system. For a direction \(\xi\), let \((\boldsymbol{\zeta}_t,\pi_t)\in {H_0^1(\Omega;\mathbb{R}^d)}\times Q\) be the solution of the linearized Stokes system \[\int_\Omega \widetilde{\mu}(\eta_t)\,\varepsilon(\boldsymbol{\zeta}_t):\varepsilon(\boldsymbol{v})\,\mathrm{d}x - \int_\Omega \pi_t \nabla\cdot \boldsymbol{v}\,\mathrm{d}x = - \int_\Omega \widetilde{\mu}'(\eta_t) \varepsilon({\boldsymbol{u}}[\eta_t]):\varepsilon(\boldsymbol{v})\,\xi\,\mathrm{d}x,\] \[\int_\Omega q\,\nabla\cdot\boldsymbol{\zeta}_t\,\mathrm{d}x=0 .\] Since both \({\boldsymbol{u}}[\eta_t]\) and \(\boldsymbol{\zeta}_t\) are divergence-free, the pressure terms vanish when testing with \({\boldsymbol{u}}[\eta_t]\) or \(\boldsymbol{\zeta}_t\). Therefore \[\frac{d^2}{dt^2}\widetilde{\mathcal{J}}_{\boldsymbol{u}}(\eta_t) = \omega\int_\Omega \widetilde{\mu}(\eta_t)|\varepsilon(\boldsymbol{\zeta}_t)|^2\,\mathrm{d}x - \frac{\omega}{2} \int_\Omega \widetilde{\mu}''(\eta_t)|\varepsilon({\boldsymbol{u}}[\eta_t])|^2\xi^2\,\mathrm{d}x.\] The same estimates used in the scalar case lead to the convexity of \(H\).

4 Finite element discretization↩︎

In this section, we describe our finite element discretization strategy to approximate the minimum of the auxiliary energy. We denote by \(C_h\), \(U_h\), \(\boldsymbol{U}_h\), and \(Q_h\) the finite element function spaces associated with the discretization of \({L^\Phi(\Omega)}\), \({H_0^1(\Omega)}\), \({H_0^1(\Omega;\mathbb{R}^d)}\), and \(Q\) respectively, and we assume the following

Moreover, we use the common notation \[W_h := \begin{cases} U_h, &\text{scalar case}\\ \boldsymbol{U}_{h,\mathrm{div}}:=\displaystyle\left\{{\boldsymbol{u}}_h \in \boldsymbol{U}_h:\int_\Omega\nabla\cdot{\boldsymbol{u}}\,q \,\mathrm{d}x= 0,\, \forall q \in Q_h\right\}, &\text{vector-valued case} \end{cases}\]

Instead of directly discretizing the method of minimizing movements 47 (see Remark 9), we consider the metric gradient flow approach described in [37] (see also [38]). The derivation is performed entirely at the finite-dimensional level; at the continuous level, identifying the first variation with an \(L^2\)-gradient through the mapped functional 49 would require additional regularity of the state variable. In the finite-dimensional setting, instead, derivatives are elements of \(C_h^*\), and gradients are obtained from them only after choosing the appropriate Riesz map. Unless explicitly stated, we do not need Assumption 5 to be satisfied.

4.1 Metric gradient flow↩︎

We work on the open positive cone \[C_h^+:=\left\{c_h\in C_h:\mathop{\mathrm{ess\,inf}}_{x\in\Omega} c_h(x)>0\right\},\] that avoids the need to introduce one-sided admissible directions associated with box constraints (see Theorem 10), and consider the smooth discrete energy \[\label{eq:Eh95def} {\mathcal{E}_h}(c_h):=-\omega \mathcal{J}_c(c_h)-\omega \mathcal{J}_{u_h}(c_h),\quad c_h\in C_h^+.\tag{62}\] where \(\mathcal{J}_{u_h}\) is identical to \(\mathcal{J}_u\) except that the infimum of the Dirichlet problem is sought in \(W_h\). Throughout this section, we always assume the validity of Assumption 6 or 7.

We then derive the first variation of the discrete energy for a fixed \(c_h\in C^+_h\). We consider the scalar case first. Let \(w_h\in C_h\) such that \(\delta w_h\) is an admissible direction. Recalling that \(\mu(r)=2\Phi'(r)\), the first variation of \({\mathcal{E}_h}\) in the direction \(w_h\) is \[\label{eq:first95var95temp} \begin{align} {\mathcal{DE}_h}(c_h)[w_h] &:= \left.\frac{d}{d\delta}{\mathcal{E}_h}(c_h+\delta w_h)\right|_{\delta=0}\\ &=\int_\Omega -\omega\left( \Phi'(c_h) -\frac{1}{2}\mu(c_h) -\frac{1}{2}\mu'(c_h)c_h -\frac{1}{2}\mu'(c_h)|\nabla u_h[c_h]|^2 \right)w_h +\omega \mu(c_h)\nabla u_h[c_h]\cdot \nabla\zeta_h\,\mathrm{d}x\\ &=\int_\Omega -\omega\left( -\frac{1}{2}\mu'(c_h)c_h -\frac{1}{2}\mu'(c_h)|\nabla u_h[c_h]|^2 \right)w_h +\omega \mu(c_h)\nabla u_h[c_h]\cdot \nabla\zeta_h\,\mathrm{d}x \end{align}\tag{63}\] where \(\zeta_h\) is the variation of \(u_h[c_h]\) in the direction \(w_h\), i.e. \[\zeta_h:=\left.\frac{d}{d\delta}u_h[c_h+\delta w_h]\right|_{\delta=0}.\] Differentiating 31 yields \[\label{eq:sens95scalar} \int_\Omega \mu(c_h)\nabla\zeta_h\cdot \nabla v_h\,\mathrm{d}x +\int_\Omega \mu'(c_h) \nabla u_h[c_h]\cdot \nabla v_h \,w_h \,\mathrm{d}x= 0, \quad \forall v_h\in U_h,\tag{64}\] that, when tested with \(v_h=u_h[c_h]\), gives the identity \[\label{eq:key95identity95scalar} \int_\Omega \mu'(c_h)|\nabla u_h[c_h]|^2 w_h\,\mathrm{d}x = -\int_\Omega \mu(c_h)\nabla\zeta_h \cdot \nabla u_h[c_h]\,\mathrm{d}x.\tag{65}\] Substituting 65 in 63 , we obtain \[\label{eq:first95var95scalar} {\mathcal{DE}_h}(c_h)[w_h] = \int_\Omega \frac{-\omega\mu'(c_h)}{2} \left(|\nabla u_h[c_h]|^2-c_h\right)\,w_h\,\mathrm{d}x.\tag{66}\]

We then compute the first variation of the auxiliary energy in the incompressible vector-valued case \[{\mathcal{DE}_h}(c_h)[w_h] =\int_{\Omega} -\omega\left( -\frac{1}{2}\mu'(c_h)c_h -\frac{1}{2}\mu'(c_h)\,|\varepsilon({\boldsymbol{u}}_h[c_h])|^2 \right)w_h +\omega\mu(c_h)\,\varepsilon({\boldsymbol{u}}_h[c_h])\!:\!\varepsilon(\boldsymbol{\zeta}_h)\,\mathrm{d}x\] where \[\boldsymbol{\zeta}_h:=\left.\frac{d}{d\delta}{\boldsymbol{u}}_h[c_h+\delta w_h]\right|_{\delta=0}.\] Using the notation \(\pi_h:=\left.\frac{d}{d\delta}p_h[c_h+\delta w_h]\right|_{\delta=0},\) a differentiation of the Stokes constraint 37 –38 yields \[\begin{align} \int_{\Omega}\mu(c_h)\,\varepsilon(\boldsymbol{\zeta}_h)\!:\!\varepsilon(\boldsymbol{v}_h)\,\mathrm{d}x -\int_{\Omega} \pi_h\,\nabla\cdot\boldsymbol{v}_h\,\mathrm{d}x &= -\int_{\Omega}\mu'(c_h)\,\varepsilon({\boldsymbol{u}}_h[c_h])\!:\!\varepsilon(\boldsymbol{v}_h)\,w_h\,\mathrm{d}x,\quad\forall \boldsymbol{v}_h \in \boldsymbol{U}_h,\tag{67}\\ \int_{\Omega} q_h\,\nabla\cdot\boldsymbol{\zeta}_h \,\,\mathrm{d}x&=0,\quad \forall q_h\in Q_h.\tag{68} \end{align}\] Testing the first equation with \(\boldsymbol{v}_h={\boldsymbol{u}}_h[c_h]\) and using the incompressibility of \({\boldsymbol{u}}_h[c_h]\), we arrive at \[\label{eq:first95var95vector} {\mathcal{DE}_h}(c_h)[w_h] := \int_\Omega \frac{-\omega\mu'(c_h)}{2} \left(|\varepsilon({\boldsymbol{u}}_h[c_h])|^2-c_h\right)\,w_h\,\mathrm{d}x.\tag{69}\]

The same first variation is obtained if homogeneous essential boundary conditions are replaced by homogeneous or non-homogeneous natural boundary conditions. In that case the state problem is posed on the corresponding linear space \(W_h\), normalized, as usual, when the operator has non-trivial kernel, and the right-hand side is replaced by \[F_N(v_h):=F(v_h)+\int_{\partial\Omega} g_N \,v_h\,\mathrm{d}s.\] Since \(F_N\) is independent of \(c_h\), differentiating the state equation gives the same linearized problems as above, and therefore the same cancellation. We then show how to handle the case of non-homogeneous essential boundary conditions; we just consider the scalar case, since similar reasoning holds for the incompressible vector-valued case and is omitted for brevity. With non-homogeneous essential boundary conditions, \(U_h\) is an affine space and the state problem is \[\int_\Omega \mu(c_h)\nabla u_h[c_h]\cdot \nabla v_h\,\mathrm{d}x=F(v_h), \quad \forall v_h\in U_{h,0},\] where \(U_{h,0}\) denote the space with homogeneous boundary conditions. The contribution of the elliptic problem is \[\mathcal{J}_{u_h}(c_h) = \frac{1}{2}\int_\Omega \mu(c_h)|\nabla u_h[c_h]|^2\,\mathrm{d}x-F(u_h[c_h]),\] with first variation \[\mathcal{DJ}_{u_h}(c_h)[w_h] = \frac{1}{2}\int_\Omega \mu'(c_h)|\nabla u_h[c_h]|^2w_h\,\mathrm{d}x + \int_\Omega \mu(c_h)\nabla u_h[c_h] \cdot \nabla\zeta_h\,\mathrm{d}x - F(\zeta_h).\] where \(\zeta_h\in U_{h,0}\) since the lifting is independent on \(c_h\). Therefore, \(\zeta_h\) is an admissible test function in the state equation, giving the identity \[\int_\Omega \mu(c_h)\nabla u_h[c_h]\cdot \nabla\zeta_h\,dx = \mathcal{F}(\zeta_h),\] and thus yielding the same first variation 66 .

Following [37], the differential \({\mathcal{DE}_h}(c_h)\in C_h^*\) is turned into a gradient \(\nabla_{g_{c_h}} {\mathcal{E}_h}(c_h) \in C_h\) by solving the variational problem \[g_{c_h}(\nabla_{g_{c_h}} {\mathcal{E}_h}(c_h), w_h) = {\mathcal{DE}_h}(c_h)[w_h], \quad \forall w_h \in C_h,\] where the inner product arising from the discrete metric tensor is \[\label{eq:gch95inner95product} g_{c_h}(z_h,w_h):=\bigl(\mathcal{D}\Psi(c_h)z_h,\mathcal{D}\Psi(c_h)w_h\bigr)_{{L^2(\Omega)}} = \int_\Omega\frac{\omega}{2}\mu'(c_h)z_h w_h\,\mathrm{d}x.\tag{70}\] Note that, since \(C_h\) is a finite-dimensional linear space, its tangent space at each point is canonically \(C_h\) itself. Moreover, under the sole Assumption 3, \(\omega\mu'(c_h)>0\) on \(C_h^+\), and therefore \(g_{c_h}(\cdot,\cdot)\) defines an inner product on \(C_h\) for all \(c_h \in C^+_h\). Note that this distinction between derivatives and gradients is standard in Hilbert-space PDE-constrained optimization [58] and Bayesian inverse problems, where the prior covariance determines the geometry used to identify dual quantities with admissible directions [59].

In the present setting, the metric tensor, equivalently the Riesz map \[\label{mxewzqfu} \mathcal{G}(c_h) : C_h\rightarrow C^\ast_h,\quad\langle\mathcal{G}(c_h)w_h,z_h\rangle = g_{c_h}(z_h,w_h), \quad c_h \in C^+_h,\,w_h,z_h \in C_h\tag{71}\] is explicitly known, and we do not need to solve any additional variational problem to define the gradient. This is in contrast with many gradient-flow equations, where the so-called Onsager operator \(K_h(c_h):C_h^*\to C_h\), mapping forces to velocities, is specified first and the metric tensor is only implicitly defined as its inverse \(G_h(c_h)=K_h(c_h)^{-1}\).

Using the first variation of the functional given in 66 for the scalar case, 69 for the vector-valued case, and the short-hand notation \(Du\) to denote \(\nabla u\) in the scalar case and \(\varepsilon({\boldsymbol{u}})\) in the vector-valued case, the metric gradient flow approach therefore gives \[\nabla_{g_{c_h}} {\mathcal{E}_h}(c_h) = c_h - |D u_h[c_h]|^2,\] and it allows us to consider the ordinary differential equation \[\label{eq:ode95fem} \left\{ \begin{align} \frac{dc_h}{dt} &= |D u_h[c_h]|^2 - c_h,\\ c_h(0) &= c_{h,0} \in C^+_h. \end{align} \right.\tag{72}\] Assumption 8 ensures that the pointwise ODE is equivalent to the semi-discrete problem \[\label{eq:ode95semi} \left\{ \begin{align} \int_\Omega \frac{dc_h}{dt} w_h\,\mathrm{d}x&= \int_\Omega(|D u_h[c_h]|^2 - c_h) w_h\,\mathrm{d}x, \quad \forall w_h \in C_h,\\ c_h(0) &= c_{h,0} \in C^+_h. \end{align} \right.\tag{73}\] Moreover, \({\mathcal{E}_h}\) is a Lyapunov function, since \[\frac{d}{dt}{\mathcal{E}_h}(c_h(t)) = -\int_\Omega\frac{\omega}{2}\mu'(c_h(t))\,\dot{c}_h^{\,2}(t)\,\mathrm{d}x\le 0, \quad \dot{c}_h :=\frac{dc_h}{dt}.\]

Box constraints are only needed to define the continuous functional \({\mathcal{E}_\epsilon}\). At the discrete level, however, the following result shows that they can be neglected for finite-time evolutions, and no explicit projection onto the box constraint is needed in the finite element implementation.

We prove the result in the scalar case. The vector-valued case is identical after replacing \(\nabla u_h\) by \(\varepsilon({\boldsymbol{u}}_h)\) and working in the discrete divergence-free space \(\boldsymbol{U}_{h,\mathrm{div}}\). It is omitted for brevity. We first prove that the function \[\label{eq:fode} \mathcal{F}_h:C_h^+\to C_h,\quad \mathcal{F}_h(c_h):=|\nabla u_h[c_h]|^2-c_h,\tag{74}\] is locally Lipschitz. To this end, for \(0<m<M\) we define \[K_{m,M}:=\{c_h\in C_h:\;m\le c_h\le M\},\] which is compact in \(C_h\). Moreover, for each \(c_{h,0}\in C^+_h\), the elliptic problem that defines \(\nabla u_h[c_h]\) for each \(c_h\) in the neighborhood of \(c_{h,0}\) is well-posed, since \(c_{h,0} \in K_{m,M}\) for some \(m,M\). We take two functions \(c_i\in K_{m,M}\) and use the notation \(u_i:=u_h[c_i]\), \(i=1,2\), i.e. \[\int_\Omega\mu(c_i)\nabla u_i\cdot\nabla v_h\,\mathrm{d}x=F(v_h),\quad \forall v_h\in U_h.\] Combining the equations, we obtain \[\int_\Omega\mu(c_1)\nabla(u_1-u_2)\cdot\nabla v_h\,\mathrm{d}x=\int_\Omega(\mu(c_2)-\mu(c_1))\nabla u_2\cdot\nabla v_h\,\mathrm{d}x.\] Choosing \(v_h=u_1-u_2\), and using the lower bound \[\mu(c_1)\ge \mu_{m,M}^{\min} := \min_{r\in[m,M]}\mu(r)>0,\] yields \[\mu_{m,M}^{\min}\|\nabla u_1-\nabla u_2\|_{{L^2(\Omega)}}\le\|\mu(c_2)-\mu(c_1)\|_{{L^\infty(\Omega)}}\|\nabla u_2\|_{{L^2(\Omega)}}.\] Since \(\mu\in C^1([m,M])\), \[\|\mu(c_2)-\mu(c_1)\|_{{L^\infty(\Omega)}}\le L^\mu_{m,M}\|c_2-c_1\|_{{L^\infty(\Omega)}},\quad L^\mu_{m,M}:=\max_{r\in[m,M]}|\mu'(r)|< +\infty.\] Combining with the standard estimate \[\|\nabla u_2\|_{{L^2(\Omega)}}\le \frac{\|F\|_{H^{-1}(\Omega)}}{\mu_{m,M}^{\min}},\] we obtain \[\|\nabla u_1-\nabla u_2\|_{{L^2(\Omega)}}\le L_u\|c_1-c_2\|_{{L^\infty(\Omega)}}, \quad L_u := L^\mu_{m,M}\frac{\|F\|_{H^{-1}(\Omega)}}{(\mu_{m,M}^{\min})^2},\] for all \(c_1,c_2\in K_{m,M}\). Denoting with \(C_{\mathrm{inv},h}\approx h^{-d/2}\) the constant associated with the standard inverse inequality for shape-regular meshes \(\|\nabla v_h\|_{{L^\infty(\Omega)}}\le C_{\mathrm{inv},h}\|\nabla v_h\|_{{L^2(\Omega)}},\, \forall v_h \in U_h\) (see e.g. [60]), the local Lipschitz continuity of \(\mathcal{F}_h\) then follows: \[\begin{align} \|\mathcal{F}_h(c_1)-\mathcal{F}_h(c_2)\|_{{L^\infty(\Omega)}} &\le \||\nabla u_1|^2-|\nabla u_2|^2\|_{{L^\infty(\Omega)}}+\|c_1-c_2\|_{{L^\infty(\Omega)}}\\ &\le(\|\nabla u_1\|_{{L^\infty(\Omega)}}+\|\nabla u_2\|_{{L^\infty(\Omega)}})\|\nabla u_1-\nabla u_2\|_{{L^\infty(\Omega)}}+\|c_1-c_2\|_{{L^\infty(\Omega)}}\\ &\le C_{\mathrm{inv},h}(\|\nabla u_1\|_{{L^\infty(\Omega)}}+\|\nabla u_2\|_{{L^\infty(\Omega)}})\|\nabla u_1-\nabla u_2\|_{{L^2(\Omega)}}+\|c_1-c_2\|_{{L^\infty(\Omega)}}\\ &\le \bigg(C_{\mathrm{inv},h}L_u(\|\nabla u_1\|_{{L^\infty(\Omega)}}+\|\nabla u_2\|_{{L^\infty(\Omega)}}) + 1\bigg)\|c_1-c_2\|_{{L^\infty(\Omega)}},\\ &\le L\,\|c_1-c_2\|_{{L^\infty(\Omega)}}, \end{align}\] where \[L:= 2 C_{\mathrm{inv},h}^2L_u\frac{\|F\|_{H^{-1}(\Omega)}}{\mu_{m,M}^{\min}} + 1.\]

Since \(C_h^+\) is open in the finite-dimensional space \(C_h\), and \(\mathcal{F}_h:C_h^+\to C_h\) is locally Lipschitz, the Picard–Lindelöf theorem (see e.g. [61]) gives a unique maximal solution \[c_h\in C^1([0,T_{\max});C_h), \quad c_h(t)\in C_h^+,\] for some \(0<T_{\max}\le+\infty\). For \(t < T_{\max}\), the explicit solution of the ODE is \[\label{eq:ode95expl} c_h(t)=c_{h,0}\,e^{-t}+\int_0^t e^{-(t-s)}|\nabla u_h[c_h(s)]|^2\,ds.\tag{75}\] Since \(c_{h,0}\in C^+_h\), we have \(m_0 := \mathop{\mathrm{ess\,inf}}_{x\in\Omega}c_{h,0}(x) > 0\), \(M_0 := \mathop{\mathrm{ess\,sup}}_{x\in\Omega}c_{h,0}(x) < +\infty\), and we have the lower bound \[c_h(t)\ge m_T:=m_0\,e^{-T} > 0,\quad \forall\, 0\le t \le T< T_{\max},\] and thus we cannot have finite-time extinction since the trajectories cannot approach the boundary of \(C^+_h\) in finite time. We then estimate an upper bound on \([0,T_{\max})\), and use the notation \(y(t) := \|c_h(t)\|_{L^\infty(\Omega)}\). By coercivity of the bilinear form and the inverse estimate \[\|\nabla u_h[c_h(t)]\|_{{L^\infty(\Omega)}}^2\le\left(\frac{C_{\mathrm{inv},h}\|F\|_{H^{-1}(\Omega)}}{\mu_{\min}(t)}\right)^2,\quad \mu_{\min}(t):=\mathop{\mathrm{ess\,inf}}_{x\in\Omega}\mu(c_h(x,t)) > 0.\] In the case of \(\Phi\) convex, \(\mu\) is increasing and thus \(\mu_{\min}(t) \ge \mu(m_T)\) for \(t\le T\); from 75 we obtain \[y(t) \le K_{T,\mathrm{conv}} + (M_0 - K_{T,\mathrm{conv}})e^{-t}, \quad 0\le t \le T, \quad K_{T,\mathrm{conv}}:=\left(\frac{C_{\mathrm{inv},h}\|F\|_{H^{-1}(\Omega)}}{\mu(m_T)}\right)^2.\] In the case of concave \(\Phi\), \(\mu\) is decreasing; from Lemma 2 \[\mu(c_h(t))\ge \mu(y(t)) \ge \mu(m_T)\left(\frac{m_T}{y(t)}\right)^\frac{q}{2}, \quad q:=1-{{S}^{-}_{\phi}}\in[0,1).\] Therefore \[\|\nabla u_h[c_h(t)]\|_{{L^\infty(\Omega)}}^2\le K_{T,\mathrm{conc}}\,y(t)^q,\quad K_{T,\mathrm{conc}}:=\frac{C_{\mathrm{inv},h}^2\|F\|_{H^{-1}(\Omega)}^2}{\mu(m_T)^2m^q_T}.\] Plugging the above estimate into 75 yields \[y(t)\le M_0 + K_{T,\mathrm{conc}}\int_0^t y(s)^q\,ds.\] Using Bihari’s inequality [62] with \[k := M_0,\quad M := K_{T,\mathrm{conc}},\quad \omega(u):=u^q,\quad \Omega(u) :=\frac{u^{1-q}-u_0^{1-q}}{1-q}, \quad \Omega^{-1}(z) = \left(u_0^{1-q}+(1-q)z\right)^{\frac{1}{1-q}},\] gives \[y(t) \le\left(M_0^{1-q}+K_{T,\mathrm{conc}}(1-q)\,t\right)^\frac{1}{1-q}, \quad 0\le t \le T.\] Therefore, in both regimes, for every \(T<T_{\max}\), there exists \(B_T<+\infty\) such that \[m_T\le c_h(x,t)\le B_T,\quad 0\le t\le T,\] and we obtain \[\epsilon(T):=\min\left\{m_T,\frac{1}{B_T}\right\}.\]

We then prove global existence, and assume by contradiction that \(T_{\max}<+\infty\). Defining \(m_*:=m_0e^{-T_{\max}} > 0\), since \(m_*\le m_T\), in the convex case we have \[K_{T,\mathrm{conv}} \le K_{*,\mathrm{conv}} := \left(\frac{C_{\mathrm{inv},h}\|F\|_{H^{-1}(\Omega)}}{\mu(m_*)}\right)^2 < +\infty,\] uniformly for \(T<T_{\max}\). In the concave case, using the upper bound of statement (ii) in Lemma 2 with \(m_*\le m_T\), we have \(\mu(m_T)m_T^{q/2}\ge \mu(m_*) m_*^{q/2}\), and thus \[K_{T,\mathrm{conc}}\le K_{*,\mathrm{conc}}:=\frac{C_{\mathrm{inv},h}^2\|F\|_{H^{-1}(\Omega)}^2}{\mu(m_*)^2m^q_*} < +\infty.\] Therefore, from the previous estimates we obtain uniform lower and upper bounds \[0 < m_*\le c_h(x,t)\le B_* < +\infty,\quad 0\le t<T_{\max}.\] where \[B_* := \begin{cases} \max\left\{M_0,K_{*,\mathrm{conv}}\right\}, &\Phi \text{ convex}\\ \left(M_0^{1-q}+K_{*,\mathrm{conc}}(1-q)T_{\max}\right)^\frac{1}{1-q}, &\Phi \text{ concave } \end{cases}\] Therefore \(c_h(t)\) remains in the compact set \(K_{m_*,B_*}\subset C_h^+\). This contradicts the finite-dimensional continuation alternative for ODEs [61], which says that a maximal solution with finite maximal time must leave every compact subset of the open domain \(C_h^+\). Hence \(T_{\max}=+\infty\).

We close the Section with the proof that under appropriate conditions the solution of the ODE 72 converges to the FEM solution minimizing the energies 18 and 35 , so that approximation properties are maintained. We use results from Diening and Kreuzer [8], which generalized the quasi-norm convergence estimates initially developed for the \(p\)-Laplacian by Ebmeyer and Liu [7]. In particular, we assume that the error of the Galerkin approximation of the Euler–Lagrange equations \[\label{eq:galerkin95fem} \int_\Omega A(D\bar u_h):D v_h\,\mathrm{d}x= F(v_h), \quad \forall v_h\in W_h,\tag{76}\] is measured in terms of \[\label{eq:errornorm95fem} \|V(D\bar{u}) - V(D\bar{u}_h)\|_{L^2(\Omega)},\tag{77}\] where \(\bar{u}\) is the minimizer of problem 18 or 35 , and \[A(P):=\mu(|P|^2)P=\frac{\phi'(|P|)}{|P|}P,\quad V(P):=\sqrt{\mu(|P|^2)}\,P, \quad P\in\mathbb{R}^{Nd},\] with continuous extension \(A(0)=0\) and \(V(0)=0\), and \(N\) is \(1\) or \(d\) depending if we consider the scalar or the vector-valued case. We match the notation of [8] and use \(f\lesssim g\) when \(f,g \ge0\) and \(f\le C g\), where \(C>0\) is a constant that only depends on \({{S}^{-}_{\phi}}\), \({{S}^{+}_{\phi}}\) and the \(\Delta_2\) constants of \(\phi\) and \(\phi^*\). We write \(f\simeq g\) if both \(f\lesssim g\) and \(g\lesssim f\). We also define for \(a\ge0\), the shifted N-function \(\phi_a\) by \[\label{eq:shifted95nfunc} \frac{\phi_a'(s)}{s}:=\frac{\phi'(a+s)}{a+s},\quad s>0,\quad\phi_a(0)=0,\tag{78}\] and first state the following Lemma, which recalls results from [8].

Using the statements (i) and (ii) of the previous Lemma, we can prove the following intermediate result. The proof is given in Lemma 11.

For our convergence proof, we need the metric Hessian lower bound from [37] to hold; and this is proved under sufficient conditions in the next Theorem. Estimates for \(\lambda\) are identical to the continuous case in Theorem 7, and Remark 8 translates directly to the discrete case. The proof is given in the Appendix, see Lemma 17

Introducing the notation for the squared distance in the discrete metric \[\label{eq:Dht95def} d_h(t):=g_{c_h(t)}(\dot{c}_h(t),\dot{c}_h(t))=\frac{1}{2}\int_\Omega |\mu'(c_h(t))|(c_h(t)-|Du_h[c_h(t)]|^2)^2\,\mathrm{d}x,\tag{79}\] we are then ready to state the last Theorem of this Section. The proof is provided in Lemma 18

4.2 Choice of the finite element spaces↩︎

We first consider the case of a simplicial element \(T\subset\mathbb{R}^d\), and we use \(\mathcal{P}_{k}(T)\) to denote the space of polynomials in \(d\) variables with degree at most \(k\). Given a polynomial of degree \(k\), the choice of the finite element space for the discretization of functions in \({H_0^1(\Omega)}\) is standard \[U^{k,T}_h := \big\{ u_h \in{H_0^1(\Omega)}:u_h \in C^0(\Omega), {u_h}_{|T} \in \mathcal{P}_{k}(T)\big\}.\] Given \(u_h \in U^{k,T}_h\), we have \(\nabla{u_h}_{|T} \in \mathcal{P}_{k-1}(T)\), and thus \(|\nabla{u_h}|^2\) is a polynomial of degree \(2(k-1)\) on each \(T\); a suitable discretization space for functions in \({L^\Phi(\Omega)}\) that satisfies Assumption 8 is then \[C^{k,T}_h := \big\{ c_h \in{L^\Phi(\Omega)}:{c_h}_{|T} \in \mathcal{P}_{2k-2}(T)\big\}.\] Similar considerations hold using tensor-product cells \(K\). In such a case, we use \[U^{k,K}_h := \big\{ u_h \in{H_0^1(\Omega)}:u_h \in C^0(\Omega), {u_h}_{|K} \in \mathcal{Q}_{k}(K)\big\},\] where \(\mathcal{Q}_{k}(K)\) denotes the space of polynomials in \(d\) variables with degree at most \(k\) for each variable. Since in this case \(\nabla{u_h}_{|K} \in \mathcal{Q}_{k}(K)\), a sufficient choice for the discretization space for functions in \({L^\Phi(\Omega)}\) is \[C^{k,K}_h := \big\{ c_h \in{L^\Phi(\Omega)}:{c_h}_{|K} \in \mathcal{Q}_{2k}(K)\big\}.\]

The same matching conditions hold true when \(\nabla u_h\) is replaced by the symmetric gradient \(\varepsilon({\boldsymbol{u}}_h)\); in the case of vector-valued spaces, we will consider the Taylor-Hood pairs with \(k\ge2\) \[\begin{align} \boldsymbol{U}_h^{k,T} &:= \big\{ {\boldsymbol{u}}_h \in {H_0^1(\Omega;\mathbb{R}^d)}:{\boldsymbol{u}}_h \in C^0(\Omega;\mathbb{R}^d), {{\boldsymbol{u}}_h}_{|T} \in \mathcal{P}_{k}(T)^d\big\},\\ Q_h^{k,T} &:= \big\{ p_h \in Q :p_h \in C^0(\Omega), \,{p_h}_{|T} \in \mathcal{P}_{k-1}(T)\big\}. \end{align}\] or, in case of tensor-product cells \[\begin{align} \boldsymbol{U}_h^{k,K} &:= \big\{ {\boldsymbol{u}}_h \in {H_0^1(\Omega;\mathbb{R}^d)}:{\boldsymbol{u}}_h \in C^0(\Omega;\mathbb{R}^d), {{\boldsymbol{u}}_h}_{|K} \in \mathcal{Q}_{k}(K)^d\big\},\\ Q_h^{k,K} &:= \big\{ p_h \in Q :p_h \in C^0(\Omega),\, {p_h}_{|K} \in \mathcal{Q}_{k-1}(K)\big\}. \end{align}\]

4.3 Time discretization algorithms↩︎

We first describe the time discretization for the scalar problem. The modifications required for the incompressible vector-valued case are given at the end of the Section. Throughout, the subscript \(h,n\) denotes finite element functions at the \(n\)th iteration. Also, we drop the cell and order superscripts and always write \(U_h\), \(C_h\), \(\boldsymbol{U}_h\), \(Q_h\) for the finite element spaces, assuming matching order conditions. We do not provide convergence results for the time-advancement schemes described below; rigorous convergence proofs are left for future work.

For a given initial condition \(c_{h,0}\in C_h\), we first compute the initial state \(u_{h,0}\) by solving the weighted Poisson problem: Find \(u_{h,0} \in U_h\) such that for all \(v_h \in U_h\) \[a[c_{h,0}](u_{h,0}, v_h) = F(v_h),\] where \(a[\cdot](\cdot,\cdot)\) is given by 24 . As in our previous work on biological transportation networks [63], we use the backward Euler method to discretize 73 in time and consider the nonlinear system of equations: Given \((u_{h,n},c_{h,n})\) and a current time step \(\tau_n\), find \((u_{h,n+1},c_{h,n+1})\) such that for all \((v_h,w_h) \in U_h\times C_h\) \[\label{eq:be95mm} \begin{align} \int_\Omega\frac{c_{h,n+1} - c_{h,n}}{\tau_n}\,w_h &= \int_\Omega(|\nabla u_{h,n+1}|^2 - c_{h,n+1})\,w_h \,\mathrm{d}x,\\ a[c_{h,n+1}](u_{h,n+1}, v_h) &= F(v_h). \end{align}\tag{80}\] Note that the backward Euler ensures the strict positivity of \(c_{h}\) throughout the time iteration, since \[\label{eq:be95positive} c_{h,n+1} = \frac{\tau_n}{1+\tau_n}|\nabla u_{h,n+1}|^2 + \frac{1}{1+\tau_n}c_{h,n}.\tag{81}\] Although strict positivity of \(c_h\) is preserved throughout the time iteration, the time-step nonlinear system can be difficult to solve directly, because Newton iterates need not preserve such positivity if not at convergence.

When only the steady state is of interest, one may instead consider the split iteration \[\label{eq:be95split} \begin{align} \int_\Omega\frac{c_{h,n+1} - c_{h,n}}{\tau_n}\,w_h &= \int_\Omega(|\nabla u_{h,n}|^2 - c_{h,n+1})\,w_h \,\mathrm{d}x,\\ a[c_{h,n+1}](u_{h,n+1}, v_h) &= F(v_h), \end{align}\tag{82}\] which completely decouples the two equations and, given our choice of finite element spaces, allows us to replace the first variational equation for the time advancement of \(c_{h}\) with the exact interpolation formula \[\label{eq:be95interp} c_{h,n+1} = \frac{\tau_n}{1+\tau_n}|\nabla u_{h,n}|^2 + \frac{1}{1+\tau_n}c_{h,n}.\tag{83}\]

In numerical experiments, we found that the split-iteration scheme does not always converge to the system’s steady state and may require small time steps to avoid divergence. A robust alternative we have found, for which we do not currently have any proof of reliable convergence, borrows ideas from the pseudo-transient continuation scheme [39] and its extension to differential-algebraic equations [40]. In our pseudo iteration scheme, we first compute an intermediate value for \(c_h\), denoted by \(c_{h,n+1/2}\), by using the interpolation formula 83 , and then perform a single Newton step by linearizing 80 around \((c_{h,n+1/2},u_{h,n})\) \[\begin{bmatrix} (\tau_n^{-1}+1) M & -2 B^n\\ C^n & K^n \end{bmatrix} \begin{bmatrix} \delta c_h\\ \delta u_h \end{bmatrix} = \begin{bmatrix} 0\\ g \end{bmatrix},\] where the matrix and vector entries are \[\begin{align} M_{i,j} &= \int_\Omega \phi_i\phi_j \,\mathrm{d}x, & B^n_{i,j} &= \int_\Omega \phi_i \nabla u_{h,n}\cdot \nabla\varphi_j\,\mathrm{d}x, \\ C^n_{i,j} &= \int_\Omega \mu'(c_{h,n+1/2})\,\phi_j \nabla u_{h,n}\cdot \nabla\varphi_i\,\mathrm{d}x, & K^n_{i,j} &= \int_\Omega \mu(c_{h,n+1/2})\,\nabla\varphi_i \cdot \nabla\varphi_j\,\mathrm{d}x,\\ g_i &= F(\varphi_i) - \int_\Omega \mu(c_{h,n+1/2})\,\nabla u_{h,n} \cdot \nabla\varphi_i\,\mathrm{d}x, \end{align}\] with \(\phi_i\) and \(\varphi_i\) the basis functions for \(C_h\) and \(U_h\) respectively. Note that the first component of the residual is zero because of the choice of the partial update. Unlike the fully implicit 80 and the split 82 steps, in this case the new iterate \(c_{h,n+1}=c_{h,n+1/2} + \delta c_h\) is not guaranteed to remain strictly positive and we thus need to backtrack the direction \(\delta c_h\) (together with \(\delta u_h\)) until this condition is met. Note that, in the fully nonlinear case 80 , we always start the Newton iteration from \(c_{h,n+1/2}\).

Time advancement for the incompressible vector-valued case is analogous to the scalar case, with the difference that the weighted Poisson problem is replaced by a weighted Stokes problem, and the equations for \(c_h\) involve the symmetric gradient. We state here only the fully implicit backward Euler step: Given \(({\boldsymbol{u}}_{h,n},p_{h,n},c_{h,n})\) and a current time step \(\tau_n\), find \(({\boldsymbol{u}}_{h,n+1},p_{h,n+1},c_{h,n+1})\) such that for all \((\boldsymbol{v}_h,q_h,w_h) \in \boldsymbol{U}_h\times Q_h\times C_h\) \[\begin{align} \int_\Omega\frac{c_{h,n+1} - c_{h,n}}{\tau_n}\,w_h &= \int_\Omega(|\varepsilon({\boldsymbol{u}}_{h,n+1})|^2 - c_{h,n+1})\,w_h \,\mathrm{d}x,\\ \boldsymbol{a}[c_{h,n+1}]({\boldsymbol{u}}_{h,n+1}, \boldsymbol{v}_h) - b(\boldsymbol{v}_h,p_{h,n+1})&= F(\boldsymbol{v}_h),\\ b({\boldsymbol{u}}_{h,n+1},q_h)&=0, \end{align}\] where \(\boldsymbol{a}[\cdot](\cdot,\cdot)\) and \(b(\cdot,\cdot)\) are given in 39 . The matrices and the right-hand side defining the pseudo update are the same, except that inner products are replaced by Frobenius inner products, symmetric gradients are used in place of plain gradients, and the matrix blocks accounting for the incompressibility condition of the weighted Stokes problem are included. The velocity Schur complement obtained by eliminating the \(c_h\) variable in the pseudo iteration scheme is also positive definite.

4.4 Adaptive time step selection↩︎

Since we are only interested in the steady state of the system, in our numerical experiments we will also make use of a heuristic adaptive time step strategy. The time step is chosen adaptively using numerical residual norms \[\rho_h^c(c_h,u_h):=\|R_h^{c,*}(c_h,u_h)\|_{\ell^2}, \quad \rho^f_h(c_h,u_h):=\|\mathcal{R}_h^{f,*} (c_h,u_h)\|_{\ell^2}\] where \(\mathcal{R}_h^{c,*}(c_h,u_h)\) is the assembled residual vector corresponding to the right-hand side of the ODE, and \(\mathcal{R}_h^{f,*} (c_h,u_h)\) the assembled vector including also the constraint residual. We reject the time step if \[E(u_{h,n+1})> E(u_{h,n})+10^{-4}\text{ and } \rho_{1,n}^c>\rho_{1/2,n}^c,\] where \(E\) is the original energy 20 and \(\rho_{j,n}^c:=\rho_h^c(c_{h,n+j},u_{h,n+j})\), with the convention \(u_{h,n+1/2} = u_{h,n}\). This means that we accept an increase in the energy \(E\) if it reduces the residual norm of the right-hand side of the ODE and, vice versa, accept the step if an increase in the ODE right-hand side leads to a decrease in the original energy. In case of rejection, we restore the previous iterate and set \(\tau_{n+1}=0.1\,\tau_n\).

If the step is accepted, the next time step value is selected according to the rule \[\tau_{n+1}= \begin{cases} \tau_n,& \rho_{1,n}^c\leq \rho_{1/2,n}^c,\\ \displaystyle\tau_n \min\left\{\frac{\rho^f_{1/2,n}}{\rho^f_{1,n}},2\right\},&\rho_{1,n}^c>\rho_{1/2,n}^c, \end{cases}\] where \(\rho_{j,n}^f:=\rho_h^f(c_{h,n+j},u_{h,n+j})\). The rescaling factor used after an accepted step is thus based on two residual indicators. The first is the \(c\)-component residual, which measures how close the auxiliary variable is to the fixed-point relation \(c_h=|D u_h|^2\), and that is used to decide whether the adaptive rescaling is activated. The second indicator is the fully coupled residual as in the conventional pseudo-transient continuation, which contains both the residual of the \(c_h\)-equation and the residual of the discrete state equation. Hence, the time step is enlarged up to a conservative factor of two only if the fully coupled residual decreases, but is reduced if the fully coupled residual increases.

5 Numerical experiments↩︎

In this Section, we present numerical results using the Firedrake Python package [66]; the code used for the experiments and the generation of figures and tables is public and available at https://github.com/stefanozampini/generalized_newtonian. The initial condition of the ODE 72 is always the identity, \(c_{h,0}(x)=1\). Nonlinear problems are solved using PETSc [67], [68] via the petsc4py Python bindings [69]; the nonlinear problem 76 is solved by starting from the initial condition \(u_h[c_{h,0}]\) [70], globalized by the usual cubic backtracking line-search algorithm. Linear systems are solved using the MUMPS factorization package [71]. Throughout this section, we report results for the unscaled auxiliary energy \[\mathcal{E}(c_h) := \int_\Omega \Phi(c_h) -\frac{1}{2}\mu(c_h)c_h + \frac{1}{2}|D u_h[c_h]|^2\,\mathrm{d}x- F(u_h[c_h]),\] which is equivalent, modulo \(-\omega\) scaling, to \(\mathcal{E}_h\) given in 62 .

5.1 Validation of the auxiliary energy and \(\lambda\)-convexity estimates↩︎

5.1.1 Power-law model: scalar case↩︎

In the first set of experiments, we validate the auxiliary energy and our geodesical \(\lambda\)-convexity estimates using the power-law model. For the scalar case, we consider the bump example from [26], which uses the polynomial exact solution \(u_e(x) = (x_1^2 - 1)(x_2^2 - 1)\) on the domain \(\Omega=[-1,1]^2\). As noted in [26], \(\nabla\cdot(|\nabla u_e|^{p-2}\nabla u_e) \in {L^{p'}(\Omega)}\) if and only if \(p>\sqrt{2}\), with \(p' := \frac{p}{p-1}\) the Hölder conjugate of \(p\). Still, we can obtain a well-posed problem for all \(p > 1\) by using the datum \[F(v) = \int_\Omega |\nabla u_e|^{p-2}\nabla u_e \cdot \nabla v \,\mathrm{d}x,\] since \(|\nabla u_e|^{p-2}\nabla u_e \in C^{0,\alpha}(\Omega;\mathbb{R}^d)\) with \(\alpha=\min\{p-1,1\}\), given that \(|\nabla u_e|\) vanishes at \((\pm 1, \pm1)\) and \((0,0)\), and around these points \(|\nabla u_e| \approx |x|\), and thus the right-hand side satisfies the regularity requirements from Assumption 6.

We use a structured \(32\times 32\) grid of tensor-product cells and approximation order \(k=2\). This choice allows us to exactly represent the solution and to compare the expected convergence rate of the method of minimizing movements 48 with the observed convergence rates of the time-advancement methods introduced in Section 4, namely the backward Euler method (BE) 80 , the split iteration method (Split) 82 , and the pseudo-transient continuation method (Pseudo). From Remark 8 we expect that the rate of convergence of the fully implicit time stepping strategy will be lower-bounded by the theoretical one \[\lambda_\tau = \frac{\log{(1+\lambda \tau)}}{\tau}, \quad \lambda = \begin{cases} p-1,&p<4/3\\ \frac{p}{4},& \frac{4}{3}\le p \le 4,p \ne 2,\\ 1,&p>4. \end{cases}\] We recall here that our current theory does not cover the cases \(p<4/3\) or \(p>4\); however, from Remark 12, we expect to obtain the corresponding values for \(\lambda\) asymptotically.

Figure 1 shows representative convergence curves for the energies \(E(u_{h,n})\) (blue solid lines, see eq. 18 ) and \(\mathcal{E}(c_{h,n})\) (red) as a function of time \(t\) using the BE method with a constant time step \(\tau=1\) for the cases \(p=4/3\) (left panel) and \(p=4\) (right panel) up to final time \(T=50\). As expected, \(E(u_{h,n}) \le \mathcal{E}(c_{h,n})\) for \(p=4/3\), while the inequality is reversed for \(p=4\), see Equation 34 . The time variation of the distance of the discrete conductivity from the steady state \(d(c_{h,n},c^*)\) (green dashed lines) and the difference of the energies \(|\mathcal{E}(c_{h,n}) -E(u_{h,n})|\) (magenta dotted lines) are reported in logarithmic scale using the right \(y\)-axes, confirming convergence to the exact solution. The slopes of the experimental convergence curves in the metric are also reported, indicating perfect agreement with the predicted estimates (given in parentheses).

Convergence rates using the BE, Pseudo, and Split methods are given for different values of \(p\), ranging from \(16/15\) to \(6\), and for different constant time step values in Figure 2, running the solvers up to final time \(T=50\). The Split method fails to converge for \(\tau=1\) and \(p=5,6\), while the convergence rates of the BE and Pseudo methods are very close to each other and are always lower-bounded by the theoretical estimates and the asymptotic predictions for \(p<4/3\) and \(p>4\).

5.1.2 Power-law model: incompressible vector-valued case↩︎

For the vector-valued case, we instead consider the exact solution pair \[{\boldsymbol{u}}_e(x) = 4(1-x_1^2 - x_2^2)(x_2, -x_1), \quad p_e(x) = x^2_1+ x^2_2 - \frac{1}{|\Omega|}\int_\Omega x^2_1+ x^2_2 \,\mathrm{d}x,\] on \(\Omega=B_{1,h}(0)\), a discretization of the unit disk with four thousand linear simplices, and we use \(k=3\) for the Taylor-Hood finite element pair; essential boundary conditions are imposed by sampling \({\boldsymbol{u}}_e(x)\) on the boundary points of the linear mesh. The symmetric gradient is \(|\varepsilon({\boldsymbol{u}}_e)|=4\sqrt{2}(x^2_1+ x^2_2)\) and \(\nabla\cdot(|\varepsilon({\boldsymbol{u}}_e)|^{p-2}\varepsilon({\boldsymbol{u}}_e)) \in {L^{p'}(\Omega)}\) if and only if \(p>\frac{1+\sqrt{17}}{4}\); as before, we can identify the right-hand side as \[F(\boldsymbol{v}) = \int_\Omega |\varepsilon({\boldsymbol{u}}_e)|^{p-2}\varepsilon({\boldsymbol{u}}_e) \!:\!\varepsilon(\boldsymbol{v}) \,\mathrm{d}x-\int_\Omega p_e \nabla\cdot\boldsymbol{v}\,\mathrm{d}x,\] and find that \(|\varepsilon({\boldsymbol{u}}_e)|^{p-2}\varepsilon({\boldsymbol{u}}_e) \in C^{0,\alpha}(\Omega;\mathbb{R}^d\times \mathbb{R}^d)\) with \(\alpha=\min\{2p-2,1\}\) and thus Assumption 7 is satisfied. Representative convergence curves for the cases \(p=4/3\) and \(p=4\) are given in Figure 3, while experimental convergence rates are compared with the theoretical predictions in Figure 4. We do not discuss the results in detail as they are qualitatively similar to the scalar case. In this case as well, theoretical convergence rates are confirmed, and the Split method fails to converge for \(\tau=1\) and \(p=5,6\).

5.1.3 Carreau-Yasuda model for incompressible fluids↩︎

For the next experiment, we use the Carreau–Yasuda model [41], [42] (see also [34]), for which the viscosity law is \[\mu(r)=\mu_\infty+(\mu_0-\mu_\infty)\left(1+\rho^a r^{a/2}\right)^{\frac{n-1}{a}}, \quad a>0,\quad \rho>0,\quad \mu_0>\mu_\infty\ge0,\quad n>0,\quad n\neq 1\] where \(\mu_0\) is the low-shear viscosity, \(\mu_\infty\) is the high-shear limiting viscosity, \(\rho\) is a time-scale parameter, \(a\) controls the width of the transition between regimes, and \(n\) is the power-law index. For \(0<n<1\), the model is shear-thinning, meaning that the viscosity decreases as the shear rate increases, as for example in polymer solutions, blood, paints, and biological and industrial fluids. For \(n>1\), the model is shear-thickening; the viscosity increases with the shear rate, as for example in cornstarch and water mixtures. In the intermediate regime, the model behaves like a power-law fluid.

A direct calculation shows \[\phi(t) = \frac{t^2}{2} \left[ \mu_\infty + \Delta\mu\, {}_2F_1 \left( \frac{1-n}{a}, \frac{2}{a}; 1+\frac{2}{a}; -\rho^a t^a \right) \right],\quad \Delta\mu:=\mu_0-\mu_\infty,\] where \({}_2F_1\) is the Gauss hypergeometric function. Moreover, setting \(z:=\rho^a r^{a/2}\) and \(q:=\frac{n-1}{a}\), we find \[\mu'(r) = \frac{\Delta\mu(n-1)}{2r} z(1+z)^{q-1}, \quad \mu''(r) = \frac{\Delta\mu(n-1)}{4r^2} z(1+z)^{q-2} \bigl[(a-2)+(n-3)z\bigr],\] which gives us, using \(t^2:=r\) and 10 \[\frac{t\phi''(t)}{\phi'(t)} = 1 + \frac{2r\mu'(r)}{\mu(r)} = 1+ (n-1) \frac{ \Delta\mu\,z(1+z)^{q-1} }{ \mu_\infty+\Delta\mu(1+z)^q }.\]

Therefore, for \(n<1\), \(\Phi\) is strictly concave, while for \(n>1\), \(\Phi\) is strictly convex, and we have \[1+\frac{r\mu''(r)}{2\mu'(r)} = \frac{(a+2)+(n+1)z}{4(1+z)},\quad \lambda_e:=\inf_{r>0} \left\{ 1+\frac{r\mu''(r)}{2\mu'(r)} \right\} = \min\left\{ \frac{a+2}{4}, \frac{n+1}{4} \right\}.\]

We consider the case \(\mu_\infty=0\) and, without loss of generality, simplify the model using \(\rho=1\). From \(\mu_\infty=0\), since \(0<z/(1+z)<1\), we obtain the indices \({{S}^{-}_{\phi}}=n>0\) and \({{S}^{+}_{\phi}}=1\) in the concave case, and \({{S}^{-}_{\phi}}=1\) and \({{S}^{+}_{\phi}}=n<+\infty\) in the convex case; Assumption 5 is thus not satisfied.

The structural conditions of Theorem 11 are satisfied for \(a\le 2\), \(1/3\le n\le3\), with \(n\ne1\). Assuming that asymptotically \(|D u_h[c_h(t)]|^2\le c_h(t)\) holds (see the discussion in Remark 12), we find \[\label{eq:lambda95carreau95yasuda} \lambda = \begin{cases} \lambda_e,&a\le 2, 1/3\le n\le3,\\ 1,&a>2, n>3,\\ \min\left\{n,\lambda_e\right\},&\text{otherwise}. \end{cases}\tag{84}\] On the other hand, the requirements of Lemma 6 are not satisfied; nevertheless the numerical experiments described next show convergence of the ODE to the FEM solution, indicating that Theorem 13 can be improved. We use the same experimental setting of Section 5.1.2 in terms of manufactured solution, domain, and finite-element spaces. In this case, we use the right-hand side \[F(v) = \int_\Omega {\boldsymbol{f}}_e\cdot \boldsymbol{v}\,\mathrm{d}x, \quad {\boldsymbol{f}}_e :=-\nabla\cdot\left(\mu(|\varepsilon({\boldsymbol{u}}_e)|^2) \varepsilon({\boldsymbol{u}}_e)\right) + \nabla p_e \in L^\infty(\Omega;\mathbb{R}^d).\]

In Table 1, we report the experimental convergence rates of the BE and Pseudo method with constant time steps \(\tau=1\) measured using 79 , and compare them with the theoretical lower bound with \(\lambda\) given in 84 . In all cases, the observed rates for the BE and Pseudo methods are very close to each other, and they are larger than the theoretical lower bounds; they reflect the dependence of the rates on the model parameters and are sharpest near the limiting regimes.

5.2 Comparison against Newton’s method↩︎

In this section, we compare the convergence of Newton’s method for the solution of 76 with the Pseudo iteration scheme endowed with adaptive time step selection as described in Section 4.4. As a common stopping criterion, we use the residual of 76 , i.e., the linear form \[\mathcal{R}^*[u_h](v_h):=\int_\Omega \mu(|D u_h|^2)D u_h:D v_h\,\mathrm{d}x- F(v_h),\] and declare convergence when the Euclidean norm \(\|\mathcal{R}^*[u_h]\|_{\ell^2} < 1.\)E-8 within a maximum of 200 iterations.

We emphasize that, at present, we do not have a convergence theory for the discrete-time iteration, nor for the Pseudo method equipped with the adaptive time-stepping strategy. In the numerical experiments that follow, the initial time step \(\tau_0\) is therefore chosen empirically, with the aim of ensuring convergence while reducing the number of nonlinear iterations. Nevertheless, the results indicate that the proposed approach is robust, even for models that are not currently covered by our theoretical framework.

5.2.1 Power-law model↩︎

The first test concerns the power-law model and the scalar benchmark problem considered in [27]. We take a constant source term and impose homogeneous essential boundary conditions on the L-shaped domain \(\Omega=(-1,1)^2 \setminus [0,1)\times (-1,0]\). Table 2 reports, for several values of \(p\) ranging from \(1.5\) to \(100\), and for successive refinement levels \(r\), the number of iterations required by the Pseudo method and by Newton’s method, with the latter shown in parentheses. The initial mesh consists of approximately \(46\mathrm{K}\) triangular elements, and the polynomial degree is fixed to \(k=1\). The table also reports the initial time step \(\tau_0\); for small values of \(p\), relatively large initial time steps can be employed, whereas increasingly smaller values of \(\tau_0\) are needed to obtain a convergent iteration as \(p\) grows. The iteration counts of the two methods are comparable for \(p \le 4\). Starting from \(p=8\), however, Newton’s method fails to converge, while the Pseudo method remains convergent even in the challenging case \(p=100\). Moreover, the Pseudo method exhibits an essentially mesh-independent and \(p\)-independent number of iterations for this test case.

5.2.2 Regularized Bingham model↩︎

Here we consider the Bercovier–Engelman regularization [72] of Bingham fluids [73] for which \[\phi(t) = \nu t^2+\sigma(\sqrt{t^2+b_\epsilon^2}-b_\epsilon), \quad \mu(r) = 2\nu+\frac{\sigma}{\sqrt{r+b_\epsilon^2}},\] where \(\sigma>0\) is the so-called yield stress, \(\nu>0\) the fluid viscosity, and \(b_\epsilon>0\) a regularization parameter. The model satisfies the hypothesis of Theorem 10 in the concave regime since \[\mu'(r) =-\frac{\sigma}{2}(r+b_\varepsilon^2)^{-3/2}<0, \quad\forall r>0.\] However, it does not satisfy Assumption 5 nor the requirements of Theorem 13 since \[\frac{t\phi''(t)}{\phi'(t)} =\frac{2\nu+\sigma b_\epsilon^2(t^2+b_\epsilon^2)^{-3/2}}{ 2\nu+\sigma(t^2+b_\epsilon^2)^{-1/2}},\] and thus \({{S}^{-}_{\phi}}>0\) but \({{S}^{+}_{\phi}}=1\) (see also Remark 3). For this model, \[\mu''(r)=\frac{3\sigma}{4}(r+b_\epsilon^2)^{-5/2}, \quad 1+\frac{r\mu''(r)}{2\mu'(r)} = \frac{r+4b_\epsilon^2}{4(r+b_\epsilon^2)} .\] Hence Theorem 11 is applicable with \(\lambda=1/4\) precisely when the structural condition in statement (ii) is satisfied. Such a condition reduces to \(6\nu b_\epsilon \ge \sigma\), and it is not satisfied in the experiment we report.

We reproduce the numerical experiment described in [29], considering the minimization of the vector-valued energy in the incompressible case 35 using \(\Omega=[0,1]^2\), and the manufactured solution \[{\boldsymbol{u}}_e = (u_\infty,0), \quad u_\infty(x) := \begin{cases} 0.25\left(0.4^2 - (0.4 - 2x_2)^2\right), & 0 \le x_2 \le 0.2, \\ 0.02, & 0.2 < x_2 < 0.8, \\ 0.25\left(0.4^2 - (2x_2 - 1.6)^2\right), & 0.8 \le x_2 \le 1, \end{cases}\] imposing non-homogeneous boundary conditions to satisfy the exact solution. The experiment is performed with \(\sigma=0.3\) and \(\nu=1\), for regularization parameters \(b_\epsilon \in \{10^{-2},10^{-3},10^{-4}\}\), on a sequence of uniformly refined triangular meshes starting from an initial mesh with approximately \(8\,\mathrm{K}\) elements. We use the Taylor-Hood pair with \(k=2\) to discretize the problem in space, and the initial time step of the Pseudo method is set to \(\tau_0=1\). The iteration counts are reported in Table 3. The Pseudo method exhibits mesh-independent convergence and remains robust with respect to the regularization parameter. Moreover, it converges in a comparable, and in several cases smaller, number of iterations than Newton’s method.

5.2.3 Optimal design↩︎

In the last experiment, we consider the optimal design problem described in [74] where \[\phi(r) = \begin{cases} \mu_2 r^2/2, & 0\leq r\leq \xi_1, \\ \xi_1\mu_2\left(r-\xi_1/2\right), & \xi_1\leq r\leq \xi_2, \\ \mu_1r^2/2 - \xi_1\mu_2\left(\xi_1/2-\xi_2/2\right), & \xi_2\leq r, \end{cases}\] with parameters \[0<\xi_1<\xi_2, \quad 0<\mu_1<\mu_2, \quad \xi_1\mu_2=\xi_2\mu_1.\] Then \[\phi'(r)= \begin{cases} \mu_2r,& 0<r<\xi_1,\\ \xi_1\mu_2,& \xi_1<r<\xi_2,\\ \mu_1r,& r>\xi_2, \end{cases} ,\quad \mu(r)= \begin{cases} \mu_2,& 0<r<\xi_1^2,\\ \xi_1\mu_2 r^{-1/2},& \xi_1^2<r<\xi_2^2,\\ \mu_1,& r>\xi_2^2. \end{cases}\] \(\phi\) is only \(C^1\), and on the smooth pieces, \[\phi''(r)= \begin{cases} \mu_2, & 0<r<\xi_1,\\ 0,& \xi_1<r<\xi_2,\\ \mu_1,& r>\xi_2. \end{cases}\] Thus, in the piecewise smooth sense, \({{S}^{-}_{\phi}}=0\), \({{S}^{+}_{\phi}}=1\), and the Riesz map is degenerate because \[\mu'(r)= \begin{cases} 0, & 0<r<\xi^2_1,\\ -\frac{\xi_1\mu_2}{2}r^{-3/2},& \xi^2_1<r<\xi^2_2,\\ 0,& r>\xi^2_2. \end{cases}\] This model is an N-function, but it does not satisfy any of the other structural assumptions of our framework. We reproduce two experimental settings from [18] using the scalar model with \(k=1\), with \(\mu_1=1\), \(\mu_2=2\), and \(\xi_1 = \sqrt{2\lambda_d\mu_1/\mu_2}\) for a given \(\lambda_d>0\) and constant right-hand side; the first experiment uses \(\Omega=[0,1]^2\) and \(\lambda_d=0.0084\), while the second experiment uses an L-shaped domain \(\Omega=(-1,1)^2 \setminus [0,1)\times (-1,0]\) and \(\lambda_d = 0.0145\). Iteration counts using the Pseudo method with adaptive time-step with \(\tau_0=1\) and Newton’s method are reported in Table 4, and confirm the robustness of the Pseudo method. However, in both cases, we do not observe mesh independence in the number of iterations.

6 Conclusions↩︎

In this work we introduced an auxiliary gradient-flow framework for generalized Newtonian-type variational problems. The main idea is to replace the nonlinear constitutive dependence on the gradient or strain by an auxiliary scalar variable, so that each state solve is a uniformly elliptic weighted linear problem.

At the continuous level, we studied the auxiliary energy in a metric space adapted to the growth of the N-function \(\phi\). Under suitable assumptions, we proved lower semicontinuity, geodesical \(\lambda\)-convexity, and exponential convergence of the corresponding minimizing-movement scheme. At the discrete level, we derived a finite-dimensional metric gradient flow using an appropriate Riesz map, proved global well-posedness of the resulting ODE on the positive cone of the discretized functional space, and established convergence to the finite element solution under additional structural assumptions. At the discrete level, the choice of the Riesz map is not only a technical device for defining gradients: it is precisely what transforms the first variation of the auxiliary energy into the evolution \(\frac{dc}{dt}=|D u[c]|^2-c\).

A concrete consequence of the theory is obtained for power-law models: the auxiliary gradient-flow framework gives guaranteed convergence in the ranges \(4/3\le p < 2\) and \(2< p\le4\), while the regimes \(p<4/3\) and \(p>4\) remain outside the present theory. For these latter cases, the discrete formulation nevertheless suggests asymptotic convergence rate estimates near equilibrium, and the numerical experiments are consistent with these predictions.

We also proposed several time-discretization strategies, including an operator-splitting scheme related to the Kačanov iteration and a time-adaptive pseudo-transient method that can be implemented using scalable linear solvers. The numerical experiments show that the proposed approach is robust across several models, including power-law, Carreau–Yasuda, regularized Bingham, and optimal-design examples, and is competitive with Newton’s method in the tested regimes.

Future work will focus on refining the convergence theory of the ODE, the convergence analysis of the fully discrete algorithms, the development of adaptive mesh-refinement techniques, and possible extensions to a broader, less regular class of models, including those depending on \(\varepsilon({\boldsymbol{u}})\) [75]. The adaptive time stepping procedure can also be improved by extending the techniques in [76] and [77] to the differential-algebraic case.

Acknowledgments↩︎

SZ gratefully acknowledges Giuseppe Savaré for fruitful discussions. DB is member of the INdAM research group GNCS.

7 Appendix↩︎

Let \(g(r):=\log\phi'(r)\), \(g'(r)=\phi''(r)/\phi'(r)\); from Assumption 4 we have \({{S}^{-}_{\phi}}/r\le g'(r)\le {{S}^{+}_{\phi}}/r\). Integrating these inequalities from \(s\) to \(t\) with \(0<s\le t\) yields \[{{S}^{-}_{\phi}}(\log{t} - \log{s})\le g(t)-g(s)\le {{S}^{+}_{\phi}}(\log{t} - \log{s}).\] Exponentiating the terms \[\left(\frac{t}{s}\right)^{{{S}^{-}_{\phi}}} \le \frac{\phi'(t)}{\phi'(s)}\le \left(\frac{t}{s}\right)^{{{S}^{+}_{\phi}}},\] which is equivalent to the statement (i) for \(\phi\).

In order to prove the statement for \(\mu\), we replace \(s\) and \(t\) with their square roots and use the definition of \(\mu\) from 8 , i.e. \[\frac{\mu(s)}{\mu(t)} = \frac{\phi'(\sqrt{s})}{\phi'(\sqrt{t})}\,\frac{\sqrt{t}}{\sqrt{s}},\] then apply the bounds for \(\phi'\) with \(\sqrt{s}\le \sqrt{t}\).

Let \(r>0\) and set \(t=\sqrt r\). Then, since \(r\mu(r)=r\frac{\phi'(t)}{t}=t\phi'(t)\), it suffices to prove the statement \[\frac{1}{{{S}^{+}_{\phi}}+1}\,t\phi'(t) \le \phi(t) \le \frac{1}{{{S}^{-}_{\phi}}+1}\,t\phi'(t).\] Since \(\phi'(t) > 0\), from (i) in Lemma 2, \[\phi'(t)\left(\frac{s}{t}\right)^{{{S}^{+}_{\phi}}} \le \phi'(s)\le \phi'(t)\left(\frac{s}{t}\right)^{{{S}^{-}_{\phi}}},\quad 0<s\le t.\] Integrating from \(0\) to \(t\) yields \[\frac{1}{{{S}^{+}_{\phi}}+1}\,t\phi'(t) = \int_0^t \phi'(t)\left(\frac{s}{t}\right)^{{{S}^{+}_{\phi}}}\,\mathrm{d}s \le \phi(t)\le \int_0^t \phi'(t)\left(\frac{s}{t}\right)^{{{S}^{-}_{\phi}}}\,\mathrm{d}s =\frac{1}{{{S}^{-}_{\phi}}+1}\,t\phi'(t).\]

It remains to consider the case \(r=0\). Although \(\phi'(r)\) may only be defined for \(r>0\), the quantity \(r\phi'(r)\) admits a continuous extension at \(r=0\), since \[0 \le r\phi'(r) \le {{R}^{+}_{\phi}}\phi(r) \to 0 \quad \text{as } r\to 0^+.\] where \({{R}^{+}_{\phi}}\) is given in Assumption 2. Similar reasoning holds for \(r\mu(r)\).

We prove statement (i); statement (ii) can be obtained by combining the chain of inequalities of (i) with Lemma 8. We define for \(t>0\) \[\theta(t):=\omega\left(\frac{t\phi''(t)}{\phi'(t)}-1\right).\] From Assumption 5, we have \[\label{eq:theta95bounds} 0 < A \leq \theta(t) \leq B,\tag{85}\] with \(A = {{S}^{-}_{\phi}}-1\) and \(B = {{S}^{+}_{\phi}}-1\) when \(\Phi\) is convex, while \(A = 1-{{S}^{+}_{\phi}}\) and \(B = 1-{{S}^{-}_{\phi}}\) when \(\Phi\) is concave. From 10 we obtain \[\psi'(s)=\sqrt{\frac{\omega}{2}\mu'(s)}=\sqrt{\frac{\theta(\sqrt s)}{4}}\sqrt{\frac{\mu(s)}{s}},\] and thus, integrating from \(0\) to \(r\), and using 85 we obtain \[\label{eq:psi95temp95bounds} \frac{\sqrt{A}}{2} \int^r_0 \sqrt{\frac{\mu(s)}{s}}\,\mathrm{d}s\le \psi(r) \le \frac{\sqrt{B}}{2} \int^r_0 \sqrt{\frac{\mu(s)}{s}}\,\mathrm{d}s.\tag{86}\]

Using the upper bound for \(\mu\) from (ii) in Lemma 2 we have \[\sqrt{\frac{\mu(s)}{s}} \le \sqrt{\mu(r)}\,r^{\frac{1-{{S}^{-}_{\phi}}}{4}}\, s^{\frac{{{S}^{-}_{\phi}}-3}{4}}, \quad 0<s\le r.\] Inserting the above inequality into the upper bound in 86 , and noting that \(({{S}^{-}_{\phi}}-3)/4>-1\), yields \[\psi(r)\le \frac{\sqrt{B}}{2}\sqrt{\mu(r)}\,r^{\frac{1-{{S}^{-}_{\phi}}}{4}}\int_0^r s^{\frac{{{S}^{-}_{\phi}}-3}{4}}\,\mathrm{d}s= \frac{2\sqrt{B}}{ {{S}^{-}_{\phi}}+1}\sqrt{\mu(r)}\,r^{\frac{1-{{S}^{-}_{\phi}}}{4}} r^{\frac{{{S}^{-}_{\phi}}+1}{4}} = \frac{2\sqrt{B}}{ {{S}^{-}_{\phi}}+1}\sqrt{\mu(r)r}.\] Squaring yields \[\psi(r)^2\le C_2 r\mu(r),\quad C_2 := \begin{cases} 4\frac{{{S}^{+}_{\phi}}-1}{({{S}^{-}_{\phi}}+1)^2},&\Phi \text{ convex},\\ 4\frac{1-{{S}^{-}_{\phi}}}{({{S}^{-}_{\phi}}+1)^2},&\Phi \text{ concave}. \end{cases}\]

For the proof of the lower bound, we proceed similarly, using the lower bound for \(\mu\) in (ii) Lemma 2 \[\sqrt{\frac{\mu(s)}{s}} \ge \sqrt{\mu(r)}\,r^{\frac{1-{{S}^{+}_{\phi}}}{4}}\, s^{\frac{{{S}^{+}_{\phi}}-3}{4}}, \quad 0<s\le r,\] which yields \[\psi(r)\ge \frac{\sqrt{A}}{2}\sqrt{\mu(r)}\,r^{\frac{1-{{S}^{+}_{\phi}}}{4}}\int_0^r s^{\frac{{{S}^{+}_{\phi}}-3}{4}}\,\mathrm{d}s= \frac{2\sqrt{A}}{ {{S}^{+}_{\phi}}+1}\sqrt{\mu(r)}\,r^{\frac{1-{{S}^{+}_{\phi}}}{4}} r^{\frac{{{S}^{+}_{\phi}}+1}{4}} = \frac{2\sqrt{A}}{ {{S}^{+}_{\phi}}+1}\sqrt{\mu(r)r}.\] \[\psi(r)^2\ge C_1 r\mu(r),\quad C_1 := \begin{cases} 4\frac{{{S}^{-}_{\phi}}-1}{({{S}^{+}_{\phi}}+1)^2},&\Phi \text{ convex},\\ 4\frac{1-{{S}^{+}_{\phi}}}{({{S}^{+}_{\phi}}+1)^2},&\Phi \text{ concave}. \end{cases}\]

Let \(\{c_n\}_{n \in \mathbb{N}} \subset {L^\Phi(\Omega)}\) be a Cauchy sequence with respect to the metric \(d\). By the definition of the distance in 16 , we have \[d(c_n, c_m) = \|\Psi(c_n) - \Psi(c_m)\|_{L^2(\Omega)}.\] This implies that the sequence \(\eta_n = \Psi(c_n)\) is a Cauchy sequence in \({L^2(\Omega)}\), which is a complete metric space; thus there exists a limit function \(\eta^* \in {L^2(\Omega)}\) such that \(\eta_n \to \eta^*\) in \({L^2(\Omega)}\) and \(\eta_n\) admits a subsequence (not relabeled) for which \(\eta_n\to \eta^*\) almost everywhere in \(\Omega\) as \(n \to +\infty\).

We define the candidate limit \[c^*(x) := \psi^{-1}(\eta^*(x)) \quad \text{a.e. in } \Omega,\] which is well defined since \(\psi: \mathbb{R} \to \mathbb{R}\) is a homeomorphism. The identity \(\Psi(c^*) = \eta^* \in {L^2(\Omega)}\) ensures that \(c^* \in {L^\Phi(\Omega)}\) according to Lemma 9. Finally, by the definition of the metric, we have \[\lim_{n \to +\infty} d(c_n, c^*) = \lim_{n \to +\infty} \|\Psi(c_n) - \Psi(c^*)\|_{{L^2(\Omega)}} = \lim_{n \to +\infty} \|\eta_n - \eta^*\|_{{L^2(\Omega)}} = 0,\] and thus, the Cauchy sequence \(\{c_n\}\) converges to \(c^* \in {L^\Phi(\Omega)}\), and the space is complete.

To prove that \(\gamma_g(s)\) defined in 17 is in fact a geodesic, a simple substitution shows \[\begin{align} d(\gamma_g(s),\gamma_g(r))^2 &= \int_\Omega |\Psi(\gamma_g(s)) - \Psi(\gamma_g(r))|^2\,\mathrm{d}x\\ &= \int_\Omega |(1-s)\Psi(c_0) + s\,\Psi(c_1) - (1-r)\Psi(c_0) - r\,\Psi(c_1)|^2\,\mathrm{d}x\\ &= \int_\Omega |(r-s)\Psi(c_0) -(r-s)\Psi(c_1)|^2\,\mathrm{d}x\\ &= |r-s|^2\int_\Omega |\Psi(c_0) -\Psi(c_1)|^2\,\mathrm{d}x. \\ \end{align}\]

If \(P=0\), then the left-hand side is zero. Assume \(P\neq0\), and use the notation \(a:=|P|\), \(b:=\sqrt r\), \(z:=|a-b|\), and \(m(t):=\mu(t^2)=\phi'(t)/t\). Then \[|\mu(r)-\mu(|P|^2)|\,|P|=|m(b)-m(a)|\,a, \quad |\mu'(r)|\,(r-|P|^2)^2=|\mu'(b^2)|\,(b^2-a^2)^2,\] and the desired estimate is equivalent to \[\phi_a^*(|m(b)-m(a)|\,a)\lesssim |\mu'(b^2)|\,(b^2-a^2)^2.\] Let \[s(t):=\frac{t\phi''(t)}{\phi'(t)}, \quad {{S}^{-}_{\phi}}\le s(t)\le {{S}^{+}_{\phi}}.\] Since \[m'(t)=\frac{t\phi''(t)-\phi'(t)}{t^2}, \quad \mu'(t^2)=\frac{m'(t)}{2t},\] we have \[({{S}^{-}_{\phi}}-1)\frac{\phi'(t)}{t} \le t m'(t) \le ({{S}^{+}_{\phi}}-1)\frac{\phi'(t)}{t}, \quad \frac{{{S}^{-}_{\phi}}-1}{2}\frac{\phi'(t)}{t^3} \le \mu'(t^2) \le \frac{{{S}^{+}_{\phi}}-1}{2}\frac{\phi'(t)}{t^3}.\] The extra assumptions \({{S}^{-}_{\phi}}>1\) in the convex case and \({{S}^{+}_{\phi}}<1\) in the concave case ensure that \(s(t) -1\) is uniformly separated from zero and allow us to consider the equivalences \[\begin{align} |t m'(t)|&\simeq\frac{\phi'(t)}{t} ,\tag{87}\\ |\mu'(t^2)|&\simeq\frac{\phi'(t)}{t^3}.\tag{88} \end{align}\]

We split the proof into two cases. First we assume \(b\ge a/2\) and prove that \[\label{eq:temp95m95bound951} |m(b)-m(a)|\,a\lesssim \phi_a'(z).\tag{89}\] When \(a/2\le b\le2a\), \(z = |a-b|\le a\), and \(\rho \simeq a \simeq a+z\) for all \(\min\{a,b\}\le \rho \le \max\{a,b\}\). By the mean value theorem, using 87 \[|m(b)-m(a)|\,a\le a\int_{\min\{a,b\}}^{\max\{a,b\}} |m'(\rho)|\,d\rho \simeq \int_{\min\{a,b\}}^{\max\{a,b\}} a\frac{\phi'(\rho)}{\rho^2}\,d\rho \simeq z\frac{\phi'(a+z)}{a+z} =\phi_a'(z).\] When \(b>2a\), then \(z = b-a\) and \(z\simeq b\). In the case of convex \(\Phi\), \(m\) is increasing, and thus, using the definition of \(\phi_a'(z)\) in 78 \[|m(b)-m(a)|\,a = (m(b)-m(a))\,a\le m(b)a \le m(b)b =\phi'(b) = \frac{b}{b-a}\phi_a'(z) \le 2 \phi_a'(z) \lesssim \phi_a'(z).\] When \(\Phi\) is concave, \(m\) is decreasing, and \[|m(b)-m(a)|\,a = (m(a)-m(b))\,a\le m(a)a = \phi'(a) \le \phi'(b) \simeq \phi_a'(z).\] Thus, in both cases, 89 holds and, since \(\phi_a^*\) is increasing, the uniform \(\Delta_2\)-condition for \(\phi_a^*\) gives \[\phi_a^*(|m(b)-m(a)|a)\lesssim \phi_a^*(\phi_a'(z)) \simeq \phi_a(z),\] where we have used the standard Orlicz identity \(\phi_a^*(\phi_a'(z))\simeq \phi_a(z)\), see e.g. [8]. By (ii) in Lemma 5, \[\phi_a(z)\simeq z^2\frac{\phi'(a+z)}{a+z}.\] Therefore, since \(b\ge a/2\), we have \(a+z\simeq b\), \(a+b\simeq b\); using 88 yields \[\phi_a(z) \simeq z^2\frac{\phi'(b)}{b} \simeq|\mu'(b^2)|z^2b^2\simeq|\mu'(b^2)|\,(b^2-a^2)^2,\] which proves the desired estimate for \(b\ge a/2\).

Now assume \(b<a/2\). This is where we use the additional restrictions \({{S}^{-}_{\phi}}\ge1/3\) in the concave case and \({{S}^{+}_{\phi}}\le 3\) in the convex case. First consider the convex case; then \(m(t)\) is increasing and \[|m(b)-m(a)|\,a = (m(a)-m(b))\,a\le m(a)a =\phi'(a).\] Moreover, since by definition, we have \(\phi'_a(a) = \phi'(2a)/2\), from the growth estimate in Lemma 2 we obtain \[\phi'(a) \leq 2^{-{{S}^{-}_{\phi}}}\phi'(2a) = 2^{1-{{S}^{-}_{\phi}}}\phi_a'(a).\] Using the monotonicity of \(\phi^*_a\), \[\phi^*_a(|m(b)-m(a)|\,a) \lesssim \phi_a^*(\phi'_a(a)) \lesssim \phi_a(a) \simeq a\phi'(a),\] where the last step follows by combining (i) and (ii) in Lemma 5. Using the growth estimate from Lemma 2 again and \({{S}^{+}_{\phi}}\le3\) yields \[a\phi'(a)\le a\phi'(b)\left(\frac{a}{b}\right)^{{S}^{+}_{\phi}}\leq \frac{\phi'(b)}{b^3}a^4.\] Since \(b<a/2\) we have \((b^2-a^2)^2\simeq a^4\), and thus, using 88 \[\phi_a^*(|m(b)-m(a)|\,a) \lesssim |\mu'(b^2)|\,(b^2-a^2)^2.\] Finally, consider the concave case; since \(m(t)\) is decreasing we have \[|m(b)-m(a)|a \le m(b)a = \frac{\phi'(b)}{b}a = \lambda\phi'(b), \quad \lambda:=\frac{a}{b}>2.\] We then prove that \[(\phi_a')^{-1}(\lambda\phi'(b)) \lesssim b\lambda^{1/{{S}^{-}_{\phi}}}.\] To this end, we set \[y:=C b\lambda^{1/{{S}^{-}_{\phi}}}, \quad C:=2^{1/{{S}^{-}_{\phi}}}.\] Since \({{S}^{-}_{\phi}}\le1\) and \(\lambda>1\), we have \(y\ge C b\lambda \ge a\) and \(a+y\le 2y\), yielding \[\phi_a'(y) = y\frac{\phi'(a+y)}{a+y} \ge \frac{1}{2} \phi'(y).\] Since \(y\ge 2 b\), the two-point growth estimate for \(\phi'\) in Lemma 2 \[\phi'(y) \ge \left(\frac{y}{b}\right)^{{S}^{-}_{\phi}}\phi'(b) = C^{{S}^{-}_{\phi}}\lambda \phi'(b) = 2\lambda\phi'(b).\] Thus \[\phi_a'(y)\ge \lambda\phi'(b),\] and since \(\phi_a'\) is increasing, \[(\phi_a')^{-1}(\lambda\phi'(b)) \le y = C b\lambda^{1/{{S}^{-}_{\phi}}}.\] Using the basic inequality from 4 \[\phi_a^*(t)\le (\phi_a')^{-1}(t)t,\quad t\ge0,\] with \(t=\lambda\phi'(b)\), we get, using the extra assumption \({{S}^{-}_{\phi}}\ge1/3\) \[\phi_a^*(\lambda\phi'(b)) \lesssim b\phi'(b)\lambda^{1+1/{{S}^{-}_{\phi}}} \le b\phi'(b)\lambda^4 = \frac{\phi'(b)}{b^3}a^4.\] Since \(b<a/2\), we have \[(b^2-a^2)^2\simeq a^4.\] Using 88 we obtain \[\phi_a^*(|m(b)-m(a)|a) \lesssim |\mu'(b^2)|\,(b^2-a^2)^2.\] This proves the estimate in the concave regime.

Note that the regularity assumptions of [37] hold. We derive the estimate without distinguishing between the scalar and the vector-valued case, and, with abuse of notation, we will always use \(Du:Dv\) to denote the (Frobenius) inner product between (symmetric) gradients. Differentiating \(\mathcal{R}_h(c_h)=c_h-|Du_h[c_h]|^2\) in the direction \(z_h\in C_h\) we get \[\mathcal{DR}_h(c_h)\,z_h = \left.\frac{d}{d\delta}\mathcal{R}_h(c_h+\delta z_h)\right|_{\delta=0}= z_h - 2\,Du_h[c_h]:D\zeta_h,\quad \zeta_h:=\left.\frac{d}{d\delta}u_h[c_h+\delta z_h]\right|_{\delta=0},\] and hence \[\begin{align} \langle \mathcal{G}(c_h)z_h,\mathcal{DR}_h(c)z_h\rangle &=\int_\Omega \frac{\omega}{2}\mu'(c_h)\,z_h\,(z_h-2\,Du_h[c_h]:D\zeta_h)\,\mathrm{d}x\notag\\ &=\int_\Omega \frac{\omega}{2}\mu'(c_h)\,z_h^2\,\mathrm{d}x -\omega\int_\Omega \mu'(c_h)\,z_h\,Du_h[c_h]:D\zeta_h\,\mathrm{d}x. \label{eq:first95part} \end{align}\tag{90}\] Furthermore \[\label{eq:second95part} \frac{1}{2}\langle \mathcal{DG}(c_h)[\mathcal{R}_h(c_h)]z_h,z_h\rangle =\int_\Omega \frac{\omega}{4}\,\mu''(c_h)\,\mathcal{R}_h(c_h)\,z_h^2\,\mathrm{d}x.\tag{91}\] Combining ?? , 90 , and 91 yields \[\label{eq:B95expand} \mathcal{B}(c_h,z_h)=\mathcal{A}(c_h,z_h) -\omega\int_\Omega \mu'(c_h)\,z_h\,Du_h[c_h]:D\zeta_h\,\mathrm{d}x +\int_\Omega \frac{\omega}{4}\,\mu''(c_h)\,\mathcal{R}_h(c_h)\,z_h^2\,\mathrm{d}x.\tag{92}\] Using 64 or 67 (depending on the scalar or vector case) tested with \(\zeta_h\), we get \[\label{eq:auc} -\int_\Omega \mu'(c_h)\,z_h\,Du_h[c_h]:D\zeta_h\,\mathrm{d}x = \int_\Omega \mu(c_h)\,|D\zeta_h|^2\,\mathrm{d}x,\tag{93}\] and arrive at \[\label{eq:B95decomp} \mathcal{B}(c_h,z_h)=\mathcal{A}(c_h,z_h) +\omega\int_\Omega \mu(c_h)\,|D\zeta_h|^2\,\mathrm{d}x +\int_\Omega \frac{\omega}{4}\,\mu''(c_h)\,\mathcal{R}_h(c_h)\,z_h^2\,\mathrm{d}x.\tag{94}\]

Now consider the convex case for \(\Phi\), for which \(\omega=1\). Then, plugging \[\omega\int \mu(c_h)|D\zeta_h|^2\,\mathrm{d}x\ge 0,\quad\omega\mu''(c_h)\,\mathcal{R}_h(c_h) = \mu''(c_h)\,(c_h-|Du[c_h]|^2),\] into 94 , yields \[\begin{align} \mathcal{B}(c_h,z_h)&\ge\int_\Omega \frac{\mu'(c_h)}{2}\left(1+\frac{c_h\,\mu''(c_h)}{2\mu'(c_h)} -\left(\frac{\mu''(c_h)}{2\mu'(c_h)}\right)_+|Du_h[c_h]|^2\right)z_h^2\,\mathrm{d}x\\ &\ge\lambda_{c_h}\int_\Omega \frac{\mu'(c_h)}{2}\,z_h^2\,\mathrm{d}x\\ &=\lambda_{c_h}\,\mathcal{A}(c_h,z_h), \end{align}\] where \[\label{eq:convex95lambdach} \lambda_{c_h} = \mathop{\mathrm{ess\,inf}}_{x\in\Omega}\left\{1+\frac{c_h(x)\,\mu''(c_h(x))}{2\mu'(c_h(x))} -\left(\frac{\mu''(c_h(x))}{2\mu'(c_h(x))}\right)_+|Du_h[c_h](x)|^2\right\}.\tag{95}\] We can then use the same arguments as in Theorem 7 and hence the thesis for (i) follows.

We then prove (ii). To this end, we need the following inequality \[\label{eq:CS95phi95prime} \int_\Omega \mu(c_h)\,|D\zeta_h|^2\,\mathrm{d}x \le \int_\Omega \frac{\mu'(c_h)^2}{\mu(c_h)}\,|Du_h[c_h]|^2\,z_h^2\,\mathrm{d}x,\tag{96}\] that can be proven with the same arguments used in 60 . Plugging 96 into 94 , using \(\omega=-1\), \[\begin{align} B_h(c_h,z_h) &\ge \int_\Omega \frac{-\mu'(c_h)}{2}\left( 1+\frac{c_h\mu''(c_h)}{2\mu'(c_h)} \right) z_h^2\,dx + \int_\Omega \left( \frac{1}{4}\mu''(c_h) - \frac{\mu'(c_h)^2}{\mu(c_h)} \right) |D u_h[c_h]|^2z_h^2\,\mathrm{d}x\\ &= \int_\Omega \frac{-\mu'(c_h)}{2}\left[ 1+\frac{c_h\mu''(c_h)}{2\mu'(c_h)} + \left( \frac{2\mu'(c_h)}{\mu(c_h)} - \frac{\mu''(c_h)}{2\mu'(c_h)} \right)|D u_h[c_h]|^2 \right] z_h^2\,\mathrm{d}x\\ &\ge \int_\Omega \frac{-\mu'(c_h)}{2}\left[ 1+\frac{c_h\mu''(c_h)}{2\mu'(c_h)} - \left( \frac{\mu''(c_h)}{2\mu'(c_h)} - \frac{2\mu'(c_h)}{\mu(c_h)} \right)_+|D u_h[c_h]|^2 \right] z_h^2\,\mathrm{d}x\\ &\ge \lambda_{c_h}\mathcal{A}(c_h,z_h), \end{align}\] where \[\label{eq:concave95lambdach} \lambda_{c_h} = \mathop{\mathrm{ess\,inf}}_{x\in\Omega}\left\{ 1+\frac{c_h(x)\mu''(c_h(x))}{2\mu'(c_h(x))} - \left( \frac{\mu''(c_h(x))\mu(c_h(x))-4\mu'(c_h(x))^2}{2\mu'(c_h(x))\mu(c_h(x))} \right)_+ |D u_h[c_h](x)|^2\right\}.\tag{97}\] Again, using the same arguments as in Theorem 7, statement (ii) then follows.

Using the notation \[c_t:=c_h(t),\quad u_t:=u_h[c_h(t)],\quad\mathcal{R}_t:= c_t - |D u_t|^2.\] we have \[\dot{c}_h(t)=-\mathcal{R}_t,\quad d_h(t)=\langle \mathcal{G}(c_t)\mathcal{R}_t,\mathcal{R}_t\rangle,\] where \(d_h(t)\) is defined in 79 . Differentiating \(d_h(t)\) yields \[\begin{align} d_h'(t)&=\langle \mathcal{DG}(c_t)[\dot{c}_h(t)]\mathcal{R}_t,\mathcal{R}_t\rangle + 2\langle \mathcal{G}(c_t)\mathcal{R}_t,\mathcal{DR}_t(c_t)\dot{c}_h(t)\rangle \\ &= -\langle \mathcal{DG}(c_t)[\mathcal{R}_t]\mathcal{R}_t,\mathcal{R}_t\rangle - 2\langle \mathcal{G}(c_t)\mathcal{R}_t,\mathcal{DR}_t(c_t)\mathcal{R}_t\rangle \\ &= -2 \left( \frac{1}{2} \langle \mathcal{DG}(c_t)[\mathcal{R}_t]\mathcal{R}_t,\mathcal{R}_t\rangle + \langle \mathcal{G}(c_t)\mathcal{R}_t,\mathcal{DR}_t(c_t)\mathcal{R}_t\rangle \right) \\ &= -2\mathcal{B}(c_t,\mathcal{R}_t), \end{align}\] where \(\mathcal{B}\) is given in Equation ?? . Under the assumptions of Theorem 11, we have the metric Hessian lower bound, \[\mathcal{B}(c_t,\mathcal{R}_t)\ge \lambda g_{c_t}(\mathcal{R}_t,\mathcal{R}_t)=\lambda d_h(t),\] and thus \[d_h'(t)\le -2\lambda d_h(t).\] Therefore, by Gronwall’s inequality, \[\label{eq:temp95D95h95dec} d_h(t)\le e^{-2\lambda t}d_h(0).\tag{98}\]

Recall that \(u_t\) is the solution of the weighted linear problem \[\int_\Omega \mu(c_t)Du_t:Dv_h\,\mathrm{d}x=F(v_h),\quad \forall v_h\in W_h,\] and \(\bar{u}_h\) the solution of the nonlinear problem 76 \[\int_\Omega A(D\bar{u}_h):Dv_h\,\mathrm{d}x=F(v_h),\quad \forall v_h\in W_h.\] Choosing \(v_h=u_t-\bar u_h\) yields \[\int_\Omega A( D\bar{u}_h):(Du_t-D\bar{u}_h)\,\mathrm{d}x = \int_\Omega \mu(c_t)Du_t:(Du_t- D\bar{u}_h)\,\mathrm{d}x.\] Multiplying by \(-1\) and adding \[\int_\Omega A(Du_t):(Du_t-D\bar{u}_h)\,\mathrm{d}x\] to both sides yields \[\begin{align} \int_\Omega (A(Du_t)-A( D\bar{u}_h)):(Du_t-D\bar{u}_h)\,\mathrm{d}x &= \int_\Omega \bigl(A(Du_t)-\mu(c_t)Du_t\bigr):(Du_t- D\bar{u}_h)\,\mathrm{d}x \\ &= \int_\Omega \bigl(\mu(|Du_t|^2)-\mu(c_t)\bigr) Du_t:(Du_t-D\bar{u}_h)\,\mathrm{d}x. \end{align}\] From statement (iii) in Lemma 5 \[(A(Du_t)-A(D\bar{u}_h)):(Du_t-D\bar{u}_h) \simeq |V(Du_t)-V(D\bar{u}_h)|^2,\] and therefore there exists \(c_0>0\) such that \[c_0 \|V(Du_t)-V(D\bar{u}_h)\|_{L^2(\Omega)}^2 \le \int_\Omega |\mu(|Du_t|^2)-\mu(c_t)| |Du_t| |Du_t-D\bar{u}_h|\,\mathrm{d}x.\] From the shifted Young inequality in statement (iv) of Lemma 5 with shift \(a:=|Du_t|\), for every \(\delta>0\), \[|\mu(|Du_t|^2)-\mu(c_t)| |Du_t| |Du_t-D\bar{u}_h| \le \delta \phi_{|Du_t|}(|Du_t-D\bar{u}_h|) + C_\delta \phi^*_{|Du_t|} \left( |\mu(|Du_t|^2)-\mu(c_t)| |Du_t| \right),\] and thus \[\begin{align} c_0 \|V(Du_t)-V(D\bar{u}_h)\|_{L^2(\Omega)}^2 &\le \delta \int_\Omega \phi_{|Du_t|}(|Du_t-D\bar{u}_h|)\,\mathrm{d}x + C_\delta \mathcal{C}_h(t)\\ &\le C_0\,\delta \|V(Du_t)-V(D\bar{u}_h)\|_{L^2(\Omega)}^2 + C_\delta \mathcal{C}_h(t), \end{align}\] where \[\mathcal{C}_h(t) := \int_\Omega \phi^*_{|D u_t|} \left( |\mu(c_t)-\mu(|Du_t|^2)| |D u_t| \right)\,\mathrm{d}x,\] and the last step follows from the shifted natural-distance equivalence in statement (iii) of Lemma 5. Choosing \(\delta\) sufficiently small, for example \(C_0 \delta \le c_0/2\), we obtain \[\label{eq:temp95VDU} \|V(Du_t)-V(D\bar{u}_h)\|_{L^2(\Omega)}^2 \le \frac{2C_\delta}{c_0} \mathcal{C}_h(t).\tag{99}\] By Lemma 6 we have the pointwise estimate \[\phi^*_{|Du_t|} \left( |\mu(c_t)-\mu(|Du_t|^2)|\,|Du_t| \right) \lesssim |\mu'(c_t)|\,|c_t-|Du_t|^2|^2.\] Hence, using 98 \[\mathcal{C}_h(t) \lesssim \int_\Omega |\mu'(c_t)| (c_t-|Du_t|^2)^2\,\mathrm{d}x = 2d_h(t)\lesssim e^{-2\lambda t}d_h(0).\] Combining this with 99 and taking the square root, we conclude the proof.

For compactness, with the notation from Section 4.3, set \(\widehat c_h:=c_{h,n+1/2}\), \(\widehat u_h:=u_{h,n}\), and \(\alpha_n:=\frac{\tau_n}{1+\tau_n}\). We next identify the bilinear form represented by \(S^n\) given in ?? . Let \(P_h:{L^2(\Omega)}\to C_h\) denote the \(L^2\)-projection onto \(C_h\), defined by \[\int_\Omega P_h f\,w_h\,\mathrm{d}x=\int_\Omega f w_h\,\mathrm{d}x,\quad \forall w_h\in C_h.\] For \(z_h\in U_h\), the definition of \(B^n\) implies \[M^{-1}B^n z_h= P_h(\nabla\widehat u_h\cdot \nabla z_h).\] Therefore, for \(z_h,v_h\in U_h\), \[\langle S^n z_h,v_h\rangle= \int_\Omega \mu(\widehat c_h)\nabla z_h\cdot\nabla v_h\,\mathrm{d}x + 2\alpha_n\int_\Omega\mu'(\widehat c_h)P_h(\nabla\widehat u_h\cdot\nabla z_h)(\nabla\widehat u_h\cdot\nabla v_h)\,\mathrm{d}x.\] By Assumption 8 \[\nabla z_h\cdot\nabla v_h= \frac{1}{2} \left( |\nabla(z_h+v_h)|^2-|\nabla z_h|^2-|\nabla v_h|^2\right)\in C_h\] for all \(z_h,v_h\in U_h\). Therefore the \(L^2\)-projection is exact \(P_h(\nabla\widehat u_h\cdot\nabla z_h)=\nabla\widehat u_h\cdot\nabla z_h\), and the Schur complement is represented by the symmetric bilinear form \[s^n(z_h,v_h)=\int_\Omega\mu(\widehat c_h)\nabla z_h\cdot\nabla v_h\,\mathrm{d}x + 2\alpha_n\int_\Omega\mu'(\widehat c_h)(\nabla\widehat u_h\cdot\nabla z_h) (\nabla\widehat u_h\cdot\nabla v_h)\,\mathrm{d}x.\]

We now prove the positive definiteness and take \(z_h=v_h\) \[s^n(v_h,v_h)=\int_\Omega\mu(\widehat c_h)|\nabla v_h|^2\,\mathrm{d}x + 2\alpha_n \int_\Omega\mu'(\widehat c_h)(\nabla\widehat u_h\cdot\nabla v_h)^2\,\mathrm{d}x.\] When \(\Phi\) is convex, \(\mu'(\widehat c_h)\geq 0\), and the second term is nonnegative, therefore, the statement is proved since \(\mu(\widehat c_h)> 0\). When \(\Phi\) is concave \(\mu'(\widehat c_h)<0\); from 83 we have the pointwise inequality \(\alpha_n|\nabla\widehat u_h|^2 \leq \widehat c_h.\) Therefore, since by Assumption 4 \[\mu(r)+2r\mu'(r)= \frac{\sqrt r\,\phi''(\sqrt r)}{\phi'(\sqrt r)} \mu(r)\ge {{S}^{-}_{\phi}}\mu(r)>0,\quad\forall r > 0\] we obtain \[\begin{align} \mu(\widehat c_h)|\nabla v_h|^2+2\alpha_n\mu'(\widehat c_h)(\nabla\widehat u_h\cdot\nabla v_h)^2 &\ge \mu(\widehat c_h)|\nabla v_h|^2+2\alpha_n\mu'(\widehat c_h)|\nabla\widehat u_h|^2|\nabla v_h|^2\\ &\ge \left( \mu(\widehat c_h)+2\widehat c_h\mu'(\widehat c_h)\right)|\nabla v_h|^2\\ &\ge {{S}^{-}_{\phi}}\mu(\widehat c_h)|\nabla v_h|^2, \end{align}\] and thus \[s^n(v_h,v_h)\geq {{S}^{-}_{\phi}}\int_\Omega\mu(\widehat c_h)|\nabla v_h|^2\,\mathrm{d}x>0,\] and the statement is proved.

References↩︎