3.1. Droplet-dynamics-induced distances and MFCs

The energy dissipation law (1.2) and its associated Fisher information functional (1.3) naturally induce a metric distance between two positive surface heights $h_0$ and $h_1$, as we define in the following.

Definition 3.1 (Distance functional)

Define a distance functional $\mathrm {Dist}_{V_1,V_2}\colon \mathcal {M}\times \mathcal {M}\rightarrow \mathbb {R}_+$ as below, where the space $\mathcal {M} = \{h\in L^1(\varOmega ): h\geq 0\}$. Consider the following optimal control problem:

(3.1a)\begin{equation} \mathrm{Dist}_{V_1, V_2}(h_0, h_1)^2:=\inf_{h, \boldsymbol{v}_1, v_2}\quad\int_0^1\int_\varOmega \left(|\boldsymbol{v}_1|^2 V_1(h)+ |v_2|^2V_2(h)\right)\,{\rm d}\kern0.7pt x\,{\rm d}t, \end{equation}

where the infimum is taken among $h(t,x)\colon [0,1]\times \varOmega \rightarrow \mathbb {R}_+$, $\boldsymbol {v}_1(t,x)\colon {[0,1]}\times \varOmega \rightarrow \mathbb {R}^d$, $v_2(t,x)\colon [0,1]\times \varOmega \rightarrow \mathbb {R}$, such that $h$ satisfies a reaction–diffusion-type equation with drift vector field $\boldsymbol {v}_1$, drift mobility $V_1$, reaction rate $v_2$, reaction mobility $V_2$, connecting initial and terminal surface heights $h_0$, ${h_1}\in \mathcal {M}$:

(3.1b)\begin{equation} \left. \begin{gathered} \partial_t h + \boldsymbol{\nabla}\boldsymbol{\cdot}( V_1(h) \boldsymbol{v}_1)=V_2(h)v_2,\quad (t,x)\in {[0,1]}\times \varOmega,\\ h(0, x)=h_0(x),\quad {h(1,x)=h_1(x)}, \end{gathered} \right\} \end{equation}

with no-flux boundary condition $V_1(h)\boldsymbol {v}_1\boldsymbol {\cdot }\boldsymbol {\nu }|_{\partial \varOmega }=0$.

Remark 3.2 (On convexity)

We illustrate here that the optimisation problem in Definition 3.1 is a non-convex optimisation problem with a linear constraint. To do so, we introduce variables $\boldsymbol {m} = V_1(h)\boldsymbol {v}_1$ and $s= V_2(h)v_2$, which make the constraint (3.1b) a linear equation: $\partial _t h + \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {m} -s = 0$. Then the objective functional in (3.1a) is

(3.2)\begin{equation} \mathcal{L}(h, \boldsymbol{m}, s) = \int_0^1\int_\varOmega \left(\frac{|\boldsymbol{m}|^2}{V_1(h)}+ \frac{|s|^2}{V_2(h)}\right)\,{\rm d}\kern0.7pt x\,{\rm d}t. \end{equation}

From the convexity of mobility functions $V_1(h)$ and $V_2(h)$, we can show that both terms ${|\boldsymbol {m}|^2}/{V_1(h)}=\sum _{i=1}^d({|m_i|^2}/{V_1(h)})$ and ${s^2}/{V_1(h)}$ fail to be convex in the respective variables. Here the vector $\boldsymbol {m} = (m_1,\ldots, m_d)$. In fact, the determinant of the Hessian matrix for the term ${|m_1|^2}/{V_1(h)}$ for the variables $(h,m_1)$ is

(3.3)\begin{equation} \det \begin{bmatrix} \left( 2\dfrac{V_1'(h)^2}{V_1(h)^{3}}- \dfrac{V_1''(h)}{V_1(h)^{2}}\right)|m_1|^2 & -2\dfrac{V_1(h)'}{V_1(h)^2}m_1 & \\ -2\dfrac{V_1(h)'}{V_1(h)^2}m_1 & 2V_1^{{-}1} \end{bmatrix} ={-}2\dfrac{V_1''(h)}{V_1(h)^{3}}|m_1|^2\le 0, \end{equation}

by positivity and convexity of $V_1(h)$. Hence, the term ${|m_1|^2}/{V_1(h)}$ is not convex.

Using the above-defined distance functional and the thin-film equation (1.1), we define the following MFC problem for droplet dynamics.

Definition 3.3 (MFC for droplet dynamics)

Given a time domain $[0,T]$, $T>0$, a potential functional $\mathcal {F}\colon \mathcal {M}\rightarrow \mathbb {R}$ and a terminal functional $\mathcal {G}\colon \mathcal {M}\rightarrow \mathbb {R}$, consider

(3.4a)\begin{equation} \inf_{h, \boldsymbol{v}_1, v_2}\int_0^T\left[\int_\varOmega \frac{1}{2}\left(|\boldsymbol{v}_1|^2 V_1(h)+ |v_2|^2V_2(h)\right)\,{{\rm d}\kern0.7pt x} - \mathcal{F}(h)\right]\,{\rm d}t+\mathcal{G}(h(T,\boldsymbol{\cdot})), \end{equation}

where the infimum is taken among $h(t,x)\colon [0,T]\times \varOmega \rightarrow \mathbb {R}_+$, $\boldsymbol {v}_1(t,x)\colon [0,T]\times \varOmega \rightarrow \mathbb {R}^d$ and $v_2(t,x)\colon [0,T]\times \varOmega \rightarrow \mathbb {R}$, such that

(3.4b)\begin{equation} \partial_t h +\boldsymbol{\nabla}\boldsymbol{\cdot}( V_1(h) \boldsymbol{v}_1)-V_2(h)v_2=\beta \left[\boldsymbol{\nabla}\boldsymbol{\cdot}(V_1(h)\boldsymbol{\nabla} P(h))-V_2(h) P(h)\right], \end{equation}

with boundary condition

(3.4c)\begin{equation} \left.V_1(h) (\boldsymbol{v}_1-\beta\boldsymbol{\nabla} P(h))\boldsymbol{\cdot}\boldsymbol{\nu}\right|_{\partial\varOmega} =0, \end{equation}

and initial surface height $h(0,\boldsymbol {\cdot }) = h_0$ in $\varOmega$. Here $\beta \ge 0$ is a non-negative number, which represents the strength of the droplet dynamics (1.1) in the constraint of MFC problem (3.4).

3.2. MFC reformulations

In this subsection, we focus on reformulations of the MFC problem in Definition 3.3 which will be suitable for a finite-element discretisation.

The first reformulation converts the constraint PDE (3.4b) to a linear constraint by a change of variables. Specifically, introducing the flux function $\boldsymbol {m}(t,x)\colon [0,T]\times \varOmega \rightarrow \mathbb {R}^d$ and source function $s(t, x)\colon [0,T]\times \varOmega \rightarrow \mathbb {R}$, such that

(3.5a,b)\begin{equation} \boldsymbol{m}=V_1(h)\left[\boldsymbol{v}_1 -\beta \boldsymbol{\nabla} P(h)\right], \quad s=V_2(h)\left[v_2-\beta P(h)\right], \end{equation}

then the MFC problem in Definition 3.3 is equivalent to the following linearly constrained optimisation problem: given a potential functional $\mathcal {F}\colon \mathcal {M}\rightarrow \mathbb {R}$ and a terminal functional $\mathcal {G}\colon \mathcal {M}\rightarrow \mathbb {R}$, consider

(3.6)\begin{align} &\inf_{h, \boldsymbol{m}, s}\int_0^T \int_\varOmega \frac12|\frac{\boldsymbol{m}}{V_1(h)}+\beta\boldsymbol{\nabla} P(h)|^2 V_1(h)\,{\rm d}\kern0.7pt x\,{\rm d}t \nonumber\\ &\qquad +\int_0^T\int_\varOmega \frac12|\frac{s}{V_2(h)}+\beta P(h)|^2V_2(h)\,{\rm d}\kern0.7pt x\,{\rm d}t -\int_0^T\mathcal{F}(h)\,{\rm d}t+\mathcal{G}(h(T,\boldsymbol{\cdot})), \end{align}

where the infimum is taken among functions $h, \boldsymbol {m}, s$, such that

(3.7)\begin{equation} \partial_t h +\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{m}-s=0 \quad\text{on } [0,T]\times \varOmega, \quad \left.\boldsymbol{m}\boldsymbol{\cdot}\boldsymbol{\nu}\right|_{[0,T]\times\partial\varOmega} =0, h(0,x) = h_0(x). \end{equation}

Expanding the product terms in (3.6), we get

(3.8)\begin{align} &\inf_{h, \boldsymbol{m}, s}\int_0^T\!\int_\varOmega \left( \frac{|\boldsymbol{m}|^2}{2V_1(h)} +\frac{|s|^2}{2V_2(h)} \right)+ \beta\left(\boldsymbol{m}\boldsymbol{\cdot}\boldsymbol{\nabla} P(h)+s\boldsymbol{\cdot} P(h) \right) \,{\rm d}\kern0.7pt x\,{\rm d}t \nonumber\\ &\qquad +\int_0^T \int_\varOmega \frac{\beta^2}2 \left(|\boldsymbol{\nabla} P(h)|^2V_1(h) +|P(h)|^2 V_2(h)\right) \,{\rm d}\kern0.7pt x\,{\rm d}t \nonumber\\ &\qquad -\int_0^T\mathcal{F}(h)\,{\rm d}t+\mathcal{G}(h(T,\boldsymbol{\cdot})). \end{align}

Using integration by parts and the constraint (3.7), we have

(3.9)\begin{align} & \int_0^T\int_\varOmega \left(\boldsymbol{m}\boldsymbol{\cdot}\boldsymbol{\nabla} P(h)+s\boldsymbol{\cdot} P(h)\right) \,{\rm d}\kern0.7pt x\,{\rm d}t \nonumber\\ &\qquad = \int_0^T \int_\varOmega P(h) \left(-\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{m}+s \right) \,{\rm d}\kern0.7pt x\,{\rm d}t \nonumber\\ &\qquad = \int_0^T\int_\varOmega P(h) \partial_t h\,{\rm d}\kern0.7pt x\,{\rm d}t= \mathcal{E}(h(T,\boldsymbol{\cdot}))- \mathcal{E}(h_0), \end{align}

where we used the definition of dynamic pressure $P(h)={\delta \mathcal {E}(h)}/{\delta h}$ in the last step. Combining these derivations and noting that $h_0$ is given, we arrive at the following equivalent formulation of the MFC problem in Definition 3.3.

Definition 3.5 (MFC reformulation I)

Consider

(3.10) \begin{align} &\inf_{h, \boldsymbol{m}, s} \int_0^T\int_\varOmega \left( \frac{|\boldsymbol{m}|^2}{2V_1(h)} +\frac{|s|^2}{2V_2(h)} \right)\,{\rm d}\kern0.7pt x\,{\rm d}t\nonumber\\ &\qquad +\int_0^T\int_\varOmega \frac{\beta^2}2 \left(|\boldsymbol{\nabla} P(h)|^2V_1(h) +|P(h)|^2 V_2(h)\right) \,{\rm d}\kern0.7pt x\,{\rm d}t \nonumber\\ &\qquad -\int_0^T\mathcal{F}(h)\,{\rm d}t+\mathcal{G}(h(T,\boldsymbol{\cdot})) +\beta \mathcal{E}(h(T,\boldsymbol{\cdot})), \end{align}

where the infimum is taken among $h$, $\boldsymbol {m}$ and $s$ satisfying (3.7).

Proposition 3.6 (MFC systems of droplet dynamics)

Let $(h, \boldsymbol {m}, s)$ be the critical point system of the MFC problem (3.10). Then there exists a function $\phi \colon [0, T]\times \varOmega \rightarrow \mathbb {R}$, such that

(3.11a,b)\begin{equation} \frac{\boldsymbol{m}(t,x)}{V_1(h(t,x))}=\boldsymbol{\nabla} \phi(t,x),\quad \frac{s(t,x)}{V_2(h(t,x))}=\phi(t,x), \end{equation}

and

(3.12) \begin{equation} \left. \begin{gathered} \partial_th(t,x) +\boldsymbol{\nabla}\boldsymbol{\cdot}(V_1(h(t,x))\boldsymbol{\nabla} \phi(t,x))-V_2(h(t,x))\phi(t,x)=0,\\ \partial_t\phi(t,x)+\frac{1}{2}\|\boldsymbol{\nabla} \phi(t,x)\|^2 V_1'(h(t,x))+\frac{1}{2}|\phi(t,x)|^2V'_2(h(t,x))\\ +\,\frac{\delta}{\delta h}\left[\mathcal{F}(h)-\frac{\beta^2}{2}\mathcal{I}(h)\right](t,x)=0, \end{gathered} \right\}\end{equation}

where $\mathcal {I}(h)$ is the generalised Fisher information functional given in (1.3), such that

(3.13) $$\begin{align} \frac{\delta}{\delta h}\mathcal{I}(h)&= \frac{1}{2}(V_1'(h)|\boldsymbol{\nabla} P(h)|^2+V_2'(h)|P(h)|^2)-\boldsymbol{\nabla}\boldsymbol{\cdot}(V_1(h)\boldsymbol{\nabla} P(h))\varPi'(h)\nonumber\\ &\quad +\alpha^2\Delta \boldsymbol{\nabla}\boldsymbol{\cdot}(V_1(h)\boldsymbol{\nabla} P(h))+\varPi'(h)V_2(h)P(h)-\alpha^2\Delta (P(h)V_2(h)), \end{align}$$

with initial and terminal time conditions

(3.14a,b)\begin{equation} h(0,x)=h_0(x),\quad \phi(T,x)={-}\frac{\delta}{\delta h}\left(\mathcal{G}(h(T,\boldsymbol{\cdot}))+\beta \mathcal{E}(h(T,\boldsymbol{\cdot}))\right). \end{equation}

We remark that a similar MFC formulation and system for reaction–diffusion equation was considered in our earlier work (Fu et al. Reference Fu, Osher and Li2023). We present the derivation of the MFC system (3.12) in Appendix A.

MFC problems and systems are generalisations of Benamou–Brenier formulae in optimal transport (Villani Reference Villani2008). This refers to setting $\beta =0$, $V_1(h)=h$ and $V_2(h)=0$. In the context of MFC of droplet dynamics, we need to address additional challenges, in which the dynamic pressure $P(h)$ given in (1.1e) involves a second-order Laplacian term. The role of this Laplacian term is well known in lubrication models, and its treatment in both forward simulations (Witelski & Bowen Reference Witelski and Bowen2003) and (boundary) optimal control contexts (Biswal et al. Reference Biswal, Ji, Elamvazhuthi and Bertozzi2024) is well understood. In addition, the treatment of the Laplacian term in controlling free interfaces has been discussed in the context of the KS equation (Tomlin & Gomes Reference Tomlin and Gomes2019) and weighted residual integral boundary layer models (Wray et al. Reference Wray, Cimpeanu and Gomes2022). However, the Laplacian term brings additional difficulties in the computation of the proposed MFC problem using finite-element methods (FEMs). This is from the fact that we need to approximate the forward–backward MFC system (3.12), in which the Laplacian term is enforced in the formulation of generalised Fisher information functional (1.3). We also comment that the dynamic pressure $P(h)$ is essential in modellingthe disjoining pressure and surface tension that govern the droplet dynamics. In numerical experiments, we demonstrate that the MFC problem with this pressure term exhibits essential patterns of droplets, including droplet spreading, transport, merging and splitting.

We introduce additional auxiliary variables to further reformulate the MFC problem (3.5). Let $\boldsymbol {n}(t,x): [0,T]\times \varOmega \rightarrow \mathbb {R}^d$, $p(t,x): [0,T]\times \varOmega \rightarrow \mathbb {R}$ and $\boldsymbol {q}(t,x): [0,T]\times \varOmega \rightarrow \mathbb {R}^d$ be defined as follows:

(3.15a–c)\begin{equation} \boldsymbol{n} = \alpha\boldsymbol{\nabla} h, \quad p ={-}\alpha\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{n}, \quad \boldsymbol{q} = \alpha \boldsymbol{\nabla} p. \end{equation}

This implies $p = -\alpha \boldsymbol {\nabla }\boldsymbol {\cdot }(\alpha \boldsymbol {\nabla } h) = -\alpha ^2\nabla ^2h$. Hence, the dynamic pressure $P(h)$ and its gradient $\boldsymbol {\nabla } P(h)$ can be expressed as follows:

(3.16a,b)\begin{equation} P(h) = U'(h)+p, \quad \boldsymbol{\nabla} P(h) = U''(h)\boldsymbol{\nabla} h + \boldsymbol{\nabla} p= \frac1{\alpha}(U''(h)\boldsymbol{n}+\boldsymbol{q}). \end{equation}

Plugging these relations back into the MFC problem (3.10), we obtain the following equivalent reformulation.

Definition 3.7 (MFC reformulation II)

Consider

(3.17)\begin{align} &\inf_{h, \boldsymbol{m}, s, \boldsymbol{n}, p, \boldsymbol{q}}\int_0^T\!\int_\varOmega \left( \frac{|\boldsymbol{m}|^2}{2V_1(h)} +\frac{|s|^2}{2V_2(h)}\right)\,{\rm d}\kern0.7pt x\,{\rm d}t\nonumber\\ &\qquad +\int_0^T\!\int_\varOmega \frac{\beta^2}{2} \left(\frac{|U''(h)\boldsymbol{n}+\boldsymbol{q}|^2}{\alpha^2}V_1(h) +|U'(h)+p|^2 V_2(h) \right) \,{\rm d}\kern0.7pt x\,{\rm d}t \nonumber\\ &\qquad -\int_0^T\mathcal{F}(h)\,{\rm d}t+\mathcal{G}(h(T,\boldsymbol{\cdot})) +\beta \int_{\varOmega}\left(U(h(T,x)) + \frac{|\boldsymbol{n}(T,x)|^2}{2}\right)\,{{\rm d}\kern0.7pt x} , \end{align}

where the infimum is taken among $h$, $\boldsymbol {m}$ and $s$ satisfying (3.7) and $\boldsymbol {n}$, $p$ and $\boldsymbol {q}$ satisfying (3.14a,b).

To simplify the notation, we collect the variables into a vector

(3.18)\begin{equation} \boldsymbol{u} := (h, \boldsymbol{m}, s, \boldsymbol{n}, p, \boldsymbol{q}), \end{equation}

and introduce $\boldsymbol {u}_T:= (h_T, \boldsymbol {n}_T)$, where $h_T:\varOmega \rightarrow \mathbb {R}_+$ is the terminal surface height and $\boldsymbol {n}_T:\varOmega \rightarrow \mathbb {R}^d$ the scaled surface height gradient at terminal time. Hence, $\boldsymbol {u}(t,x):[0,T]\times \varOmega \rightarrow \mathbb {R}^{3d+3}$ is a space–time function with $3d+3$ components and $\boldsymbol {u}_T(x):\varOmega \rightarrow \mathbb {R}^{d+1}$ is a spatial function with $d+1$ components. We further denote the functionals $H(\boldsymbol {u})$ and $H_T(\boldsymbol {u}_T)$, such that

(3.19a)\begin{align} H(\boldsymbol{u}) &:= \frac{|\boldsymbol{m}|^2}{2V_1(h)} +\frac{|s|^2}{2V_2(h)} +\frac{\beta^2}{2\alpha^2} |U''(h)\boldsymbol{n}+\boldsymbol{q}|^2V_1(h) \nonumber\\ &\quad +\frac{\beta^2}{2}|U'(h)+p|^2 V_2(h) - F(h), \end{align}

(3.19b)\begin{align} H_T(\boldsymbol{u}_T)&:= \beta\left(U(h_T) +\frac{|n_T|^2}{2}\right) +G(h_T), \end{align}

where $F(h)$ and $G(h)$ are density functions for functionals $\mathcal {F}$ and $\mathcal {G}$, i.e.

(3.20a,b)\begin{equation} \mathcal{F}(h) = \int_\varOmega F(h)\,{{\rm d}\kern0.7pt x}, \quad \mathcal{G}(h) = \int_\varOmega G(h)\,{{\rm d}\kern0.7pt x}. \end{equation}

Using this notation, the MFC problem (3.7) takes the following compact form.

Definition 3.8 (MFC problem: compact form)

Consider

(3.21a) \begin{equation} \inf_{\boldsymbol{u}, \boldsymbol{u}_T}\int_0^T\!\int_\varOmega H(\boldsymbol{u}) \,{\rm d}\kern0.7pt x\,{\rm d}t +\int_{\varOmega}H_T(\boldsymbol{u}_T)\,{{\rm d}\kern0.7pt x}, \end{equation}

subject to the constraints on the space–time domain

(3.21b)\begin{equation} \left. \begin{aligned} \partial_t h + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{m} - s & = 0,\\ -\boldsymbol{n} + \alpha \boldsymbol{\nabla} h & = 0,\\ p + \alpha\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{n} & = 0,\\ \boldsymbol{q} - \alpha\boldsymbol{\nabla} p & = 0, \end{aligned} \quad\text{on } [0,T]\times\varOmega, \right\} \end{equation}

and the constraints at terminal time

(3.21c) \begin{equation} \left. \begin{aligned} h(0,\boldsymbol{\cdot}) & = h_0,\\ h(T,\boldsymbol{\cdot}) & = h_T,\\ -\boldsymbol{n}_T + \alpha \boldsymbol{\nabla} h_T & = 0,\\ \end{aligned} \quad\text{on } \varOmega, \right\} \end{equation}

with Neumann boundary condition

(3.21d)\begin{equation} \boldsymbol{m}\boldsymbol{\cdot}\boldsymbol{\nu} = 0, \quad \text{on }[0,T]\times\partial\varOmega. \end{equation}

Remark 3.10 (On recovering the physics control variable $\zeta$)

Combining the definition in (3.5a,b) with the optimality conditions (3.11a,b) in Proposition 3.6, we get

(3.22a,b)\begin{equation} \boldsymbol{v}_1 - \beta \boldsymbol{\nabla} P(h) = \boldsymbol{\nabla} \phi, \quad v_2-\beta P(h) = \phi. \end{equation}

In particular, this implies that $\boldsymbol {v}_1$ is a gradient field, and the relation $\boldsymbol {v}_1 = \boldsymbol {\nabla } v_2$ holds. Using (2.19b) in § 2 that relate the control variables $\boldsymbol {v}_1$ and $v_2$ to the activity field $\zeta$, we obtain

(3.23)\begin{equation} \zeta = \frac{\beta P(h)+\phi}{h} = \frac{\beta(U'(h)+p)+\phi}{h}, \end{equation}

where we used the relation (3.16a,b) in the last equality.

3.3. Saddle-point problem

Finally, we reformulate the constrained optimisation problem in Definition 3.8 into a saddle-point problem using Lagrange multipliers, for which a finite-element discretisation is developed in § 4. We introduce the following four Lagrange multipliers on the space–time domain $[0,T]\times \varOmega$ for the four equations in (3.21b). They are scalar functions $\phi (t,x):[0,T]\times \varOmega \rightarrow \mathbb {R}$, $\xi (t,x):[0,T]\times \varOmega \rightarrow \mathbb {R}$, vectorial functions $\boldsymbol {\sigma }(t,x):[0,T]\times \varOmega \rightarrow \mathbb {R}^d$, $\boldsymbol {\theta }(t,x):[0,T]\times \varOmega \rightarrow \mathbb {R}^d$ and a Lagrange multiplier $\boldsymbol {\sigma }_T(x):\varOmega \rightarrow \mathbb {R}^d$ on the spatial domain (at terminal time):

(3.24)\begin{equation} \left. \begin{gathered} \partial_t h + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{m} - s = 0 \longleftrightarrow \phi,\\ -\boldsymbol{n} + \alpha \boldsymbol{\nabla} h = 0 \longleftrightarrow \boldsymbol{\sigma},\\ p + \alpha\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{n} = 0 \longleftrightarrow \xi,\\ \boldsymbol{q} - \alpha\boldsymbol{\nabla} p = 0 \longleftrightarrow \boldsymbol{\theta},\\ -\boldsymbol{n}_T + \alpha \boldsymbol{\nabla} h_T = 0 \longleftrightarrow \boldsymbol{\sigma}_T. \end{gathered} \right\} \end{equation}

Then the MFC problem (3.8) can be formulated as the following saddle-point problem:

(3.25)\begin{align} &\inf_{\boldsymbol{u}, \boldsymbol{u}_T} \sup_{\boldsymbol{\varPhi}, \boldsymbol{\sigma}_T} \int_0^T\int_\varOmega H(\boldsymbol{u}) \,{\rm d}\kern0.7pt x\,{\rm d}t +\int_{\varOmega}H_T(\boldsymbol{u}_T)\,{{\rm d}\kern0.7pt x}\nonumber\\ &\qquad +\int_0^T\int_\varOmega \left[(\partial_t h + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{m} - s)\phi +(-\boldsymbol{n} + \alpha\boldsymbol{\nabla} h)\boldsymbol{\cdot}\boldsymbol{\sigma}\right]\,{\rm d}\kern0.7pt x\,{\rm d}t\nonumber\\ &\qquad +\int_0^T\int_\varOmega \left[(p + \alpha\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{n})\xi+ (\boldsymbol{q}-\alpha\boldsymbol{\nabla} p)\boldsymbol{\cdot}\boldsymbol{\theta} \right] \,{\rm d}\kern0.7pt x\,{\rm d}t\nonumber\\ &\qquad +\int_\varOmega (-\boldsymbol{n}_T + \alpha\boldsymbol{\nabla} h_T)\boldsymbol{\cdot}\boldsymbol{\sigma}_T \,{{\rm d}\kern0.7pt x}, \end{align}

with the following boundary and initial/terminal conditions

(3.26)\begin{equation} h(0,\boldsymbol{\cdot}) = h_0, \ \ h(T,\boldsymbol{\cdot}) = h_T, \quad \text{in }\varOmega, \quad\text{and}\quad \boldsymbol{m}\boldsymbol{\cdot}\boldsymbol{\nu} = 0,\quad \text{on }[0,T]\times\partial\varOmega. \end{equation}

Here $\boldsymbol {\varPhi } = (\phi, \xi, \boldsymbol {\sigma }, \boldsymbol {\theta })$.

Next, applying integration by parts on the above saddle-point problem to move all derivatives of $\boldsymbol {u}$ to the dual variables $\boldsymbol {\varPhi }$, and using the initial and boundary conditions, we obtain

(3.27)\begin{align} &\inf_{\boldsymbol{u}, \boldsymbol{u}_T} \sup_{\boldsymbol{\varPhi}, \boldsymbol{\sigma}_T}\int_0^T\int_\varOmega H(\boldsymbol{u}) \,{\rm d}\kern0.7pt x\,{\rm d}t +\int_{\varOmega}H_T(\boldsymbol{u}_T)\,{{\rm d}\kern0.7pt x}\nonumber\\ &\qquad +\int_0^T\int_\varOmega \left[-( h \partial_t\phi+ \boldsymbol{m}\boldsymbol{\cdot}\boldsymbol{\nabla}\phi + s\phi) -(\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\sigma} + \alpha h\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma})\right]\,{\rm d}\kern0.7pt x\,{\rm d}t\nonumber\\ &\qquad +\int_0^T\int_\varOmega \left[(p \xi- \alpha\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\nabla} \xi)+ (\boldsymbol{q}\boldsymbol{\cdot}\boldsymbol{\theta}+\alpha p\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{\theta}) \right] \,{\rm d}\kern0.7pt x\,{\rm d}t\nonumber\\ &\qquad +\int_\varOmega\left[-(\boldsymbol{n}_T \boldsymbol{\cdot}\boldsymbol{\sigma}_T+ \alpha h_T \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma}_T) +h_T\phi(T, x)-h_0\phi(0, x)\right] \,{{\rm d}\kern0.7pt x}. \end{align}

In the above formulation, we assume the Lagrange multipliers $\boldsymbol {\sigma }$, $\xi$ and $\boldsymbol {\theta }$ satisfy the following Neumann boundary conditions:

(3.28)\begin{equation} \boldsymbol{\sigma}\boldsymbol{\cdot}\boldsymbol{\nu} = 0,\quad \boldsymbol{\nabla}{\xi}\boldsymbol{\cdot}\boldsymbol{\nu} = 0, \quad \boldsymbol{\theta}\boldsymbol{\cdot}\boldsymbol{\nu} = 0, \quad \text{on }[0,T]\times\partial\varOmega. \end{equation}

The saddle-point problem (3.27) is the final form of our MFC problem that is discretised in the next section. The variational structure of this problem makes the FEM an ideal candidate for such a problem. We close this section with a discussion on the proper function spaces for the primal variables $\boldsymbol {u}$ and $\boldsymbol {u}_T$ and dual variables $\boldsymbol {\varPhi }$ and $\boldsymbol {\sigma }_T$ in (3.27) which makes the integrals in (3.27) valid. The spaces are given as follows:

(3.29a) \begin{align} \boldsymbol{u}&\in \left\{ \boldsymbol{v}\in [L^2([0,T]\times \varOmega)]^{3d+3}:\quad \int_0^T\int_{\varOmega}H(\boldsymbol{v})\,{\rm d}\kern0.7pt x\,{\rm d}t <{+}\infty,\right.\nonumber\\ &\left.\qquad \text{first component of }\boldsymbol{v} \text{ is non-negative} \vphantom{\int_{\varOmega}}\right\}, \end{align}

(3.29b) \begin{align} \boldsymbol{u}_T&\in \left\{ \boldsymbol{v}_T\in [L^2(\varOmega)]^{d+1}:\quad \int_{\varOmega}H_T(\boldsymbol{v}_T)\,{{\rm d}\kern0.7pt x} <{+}\infty,\right.\nonumber\\ & \left.\qquad \text{first component of }\boldsymbol{v}_T \text{ is non-negative} \vphantom{\int_{\varOmega}}\right\}, \end{align}

(3.29c)\begin{align} \boldsymbol{\varPhi}&\in V_{\phi}\times \boldsymbol{V}_{\sigma}\times {V}_{\xi}\times \boldsymbol{V}_{\theta}, \qquad \boldsymbol{\sigma}_T\in H_0(\mathrm{div}; \varOmega), \end{align}

where

(3.29d)$$\begin{gather} V_{\phi}= H^1([0,T]\times\varOmega), \end{gather}$$

(3.29e)$$\begin{gather}\boldsymbol{V}_{\sigma}=\boldsymbol{V}_{\theta}= L^2([0,T])\otimes H_0(\mathrm{div};\varOmega), \end{gather}$$

(3.29f)$$\begin{gather}{V}_{\xi}= L^2([0,T])\otimes H_0^1(\varOmega). \end{gather}$$

Here we use the usual definition of Sobolev spaces

(3.30)$$\begin{gather} L^2(\varOmega):=\{v:\varOmega\rightarrow \mathbb{R}: \;\;\int_{\varOmega}|v|^2dx <{+}\infty\}, \end{gather}$$

(3.31)$$\begin{gather}H^1(\varOmega):=\{v\in L^2(\varOmega): \boldsymbol{\nabla} v \in [L^2(\varOmega)]^d\}, \end{gather}$$

(3.32)$$\begin{gather}H(\mathrm{div}; \varOmega):=\{\boldsymbol{v}\in [L^2(\varOmega)]^d: \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v} \in L^2(\varOmega)\}. \end{gather}$$

Moreover, $H^1_0(\varOmega )$ is the subspace of $H^1(\varOmega )$ with a zero boundary condition, and $H_0(\mathrm {div}; \varOmega )$ is the subspace of $H(\mathrm {div}; \varOmega )$ with a zero boundary condition on the normal direction.

4.1. The high-order finite-element scheme

We first partition the spatial domain $\varOmega$ into a spatial mesh $\varOmega _h=\{K_\ell \}_{\ell =1}^{N_S}$ with $N_S$ elements where each element $K_{\ell }$ is assumed to be a mapped hypercube in $\mathbb {R}^d$, and the temporal domain $[0,T]$ into a temporal mesh $I_h=\{I_j\}_{j=1}^{N_T}$ with $N_T$ segments. Denote the space–time mesh as $\varOmega _{T,h}=I_h\otimes \varOmega _h$. The function spaces in (3.29) for the saddle-point problem (3.27) indicate natural discretisation spaces for the primal and dual variables. In particular, we use the following conforming finite-element spaces to discretise the dual variables $\boldsymbol {\varPhi }$ and $\boldsymbol {\sigma }_T$:

(4.1a)$$\begin{gather} V_{\phi,h}^{k+1} = \{\psi\in V_{\phi}: v|_{I_j\times K_\ell} \in Q^{k+1}(I_j)\otimes Q^{k+1}(K_\ell)\quad \forall j, \ell\}, \end{gather}$$

(4.1b)$$\begin{gather}\boldsymbol{V}_{\sigma,h}^k = \{\boldsymbol{\tau}\in \boldsymbol{V}_{\sigma}: \boldsymbol{\tau}|_{I_j\times K_\ell} \in Q^k(I_j)\otimes RT^k(K_\ell)\quad \forall j, \ell\}, \end{gather}$$

(4.1c)$$\begin{gather}{V}_{\xi,h}^{k,k+1} = \{q\in \boldsymbol{V}_{\xi}: q|_{I_j\times K_\ell} \in Q^k(I_j)\otimes Q^{k+1}(K_\ell)\quad\forall j, \ell\}, \end{gather}$$

(4.1d)$$\begin{gather}\boldsymbol{M}_{\sigma,h}^k = \{\boldsymbol{\tau}\in H_0(\mathrm{div};\varOmega): \boldsymbol{\tau}|_{K_\ell} \in RT^k(K_\ell)\quad \forall j, \ell\}, \end{gather}$$

where $Q^k(K_\ell )$ is the tensor-product polynomial space of degree no greater than $k$ in each direction, and $RT^k(K_\ell )$ is the local Raviart–Thomas finite-element space Boffi, Brezzi & Fortin (Reference Boffi, Brezzi and Fortin2013) on the mapped hypercube $K_\ell$, for $k\ge 0$. For the primal variables $\boldsymbol {u}$ and $\boldsymbol {u}_T$, it is natural to use an integration rule space such that they are defined only on the (high-order) numerical integration points, since no derivative calculation is needed for these variables. Let $X_{t,k}:=\{\chi _{t,i}\}_{i=1}^{N_{T,q}^k}$ be the quadrature points and $\{\omega _{t,i}\}_{i=1}^{N_{T,q}^k}$ the corresponding quadrature weights on the temporal mesh $I_h$ using $(k+1)$ Gauss–Legendre (GL) integration points per line segment, and denote $X_{s,k}:=\{\chi _{s,j}\}_{j=1}^{N_{S,q}^k}$ as the quadrature points and $\{\omega _{s,j}\}_{j=1}^{N_{S,q}^k}$ as the corresponding quadrature weights on the spatial mesh $\varOmega _h$ using $(k+1)$ GL integration points per coordinate direction in each element. We approximate each component of $\boldsymbol {u}$ and $\boldsymbol {u}_T$ using the following space–time and spatial integration rule spaces, respectively:

(4.1e)\begin{equation} W_h^k :=\left\{v: X_{t,k}\times X_{s,k}\rightarrow \mathbb{R}\right\}, \quad M_h^k :=\left\{v: X_{s,k}\rightarrow \mathbb{R}\right\}. \end{equation}

Note that a function in the quadrature space $W_h^k$ can be interpreted as a vector of size $N_{T,q}^k\times N_{S,q}^k$ and a function in $M_h^k$ can be interpreted as a vector of size $N_{S,q}^k$.

Using the above finite-element spaces, we define the following discrete saddle-point problem: find the critical point of the discrete system

(4.2)\begin{align} &\inf_{\boldsymbol{u}_h, \boldsymbol{u}_{T,h}} \sup_{\boldsymbol{\varPhi}_h, \boldsymbol{\sigma}_{T,h}} \langle\!\langle H(\boldsymbol{u}_h)\rangle\!\rangle_h + \left\langle H_T(\boldsymbol{u}_{T,h})\right\rangle_h\nonumber\\ &\qquad - \langle\!\langle h_h \partial_t\phi_h+ \boldsymbol{m}_h\boldsymbol{\cdot}\boldsymbol{\nabla}\phi_h + s_h\phi_h\rangle\!\rangle_h-\langle\!\langle\boldsymbol{n}_h\boldsymbol{\cdot}\boldsymbol{\sigma}_h + \alpha h_h\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma}_h\rangle\!\rangle_h\nonumber\\ &\qquad +\langle\!\langle p_h \xi_h- \alpha\boldsymbol{n}_h\boldsymbol{\cdot}\boldsymbol{\nabla} \xi_h\rangle\!\rangle_h+ \langle\!\langle\boldsymbol{q}_h\boldsymbol{\cdot}\boldsymbol{\theta}_h+\alpha p_h\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{\theta}_h\rangle\!\rangle_h\nonumber\\ &\qquad -\left\langle\boldsymbol{n}_{T,h} \boldsymbol{\cdot}\boldsymbol{\sigma}_{T,h}+ \alpha h_{T,h} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma}_{T,h}\right\rangle_h +\left\langle h_{T,h}\phi_h(T, x)-h_0\phi_h(0, x)\right\rangle_h, \end{align}

where the variables $\boldsymbol {u}_h:=(h_h, \boldsymbol {m}_h, s_h, \boldsymbol {n}_h, p_h, q_h)\in [W_h^k]^{3d+3}$ with $h_h\ge 0$, $\boldsymbol {u}_{T,h}=(h_{T,h}, \boldsymbol {n}_{T,h})\in [M_h^{k}]^{d+1}$ with $h_{T,h}\ge 0$, $\varPhi _h:=(\phi _h,\boldsymbol {\sigma }_h, \xi _h, \boldsymbol {\theta }_h) \in V_{\phi,h}^{k+1}\times \boldsymbol {V}_{\sigma }^k \times V_{\xi, h}^{k,k+1}\times \boldsymbol {V}_{\sigma }^k$ and $\boldsymbol {\sigma }_{T,h}\in \boldsymbol {M}_{\sigma,h}^k$. Here the double bracket is the numerical integration on the space–time domain and single bracket is the numerical integration on the spatial domain defined as follows:

(4.3a,b) \begin{equation} \left\langle \,f(x)\right\rangle_h = \sum_{j=1}^{N_{s,j}^k}f(\chi_{s,j})\omega_{s,j}, \quad \langle\!\langle \,f(t, x)\rangle\!\rangle_h = \sum_{i=1}^{N_{t,i}^k}\sum_{j=1}^{N_{s,j}^k}f(\chi_{t,i}, \chi_{s,j})\omega_{t,i}\omega_{s,j}. \end{equation}

After solving for the discrete saddle-point problem (4.2), we can recover the dynamic pressure $P_h = U'(h_h)+p_h$ and the activity field $\zeta (t, \boldsymbol {x}) = ({\beta P_h + \phi _h})/{h_h}$; see Remark 3.10.

4.2. A generalised PDHG algorithm

We solve the discrete saddle-point problem (4.2) using a generalised preconditioned PDHG algorithm, which is a splitting algorithm that solve for the primal variables $\boldsymbol {u}_h$ and $\boldsymbol {u}_{T,h}$, and each component of the dual variables $\varPhi _h$ and $\boldsymbol {\sigma }_{T,h}$ sequentially. The following algorithm is a generalisation of the G-prox PDHG algorithm developed in Jacobs et al. (Reference Jacobs, Léger, Li and Osher2019).

Algorithm 1 Generalised PDHG for (4.2).

4.3. A fully discrete JKO scheme to the PDE (1.1)

Next we present a fully discrete JKO scheme and its simplified version for solving the PDE (1.1). As mentioned in Remark 3.4, by taking the functionals $\mathcal {F}=0$ and $\mathcal {G}(h) = \Delta t \mathcal {E}(h)$, and setting terminal time $T=1$ and the parameter $\beta = 0$, the MFC optimisation problem in Definition 3.3 becomes the dynamic formulation of a JKO temporal discretisation scheme which advances solution in time with step size $\Delta t$. Since the parameter $\beta =0$, we do not need the auxiliary variables $\boldsymbol {n}, p$, and $\boldsymbol {q}$ in the definition of the functional $H(\boldsymbol {u})$ in (3.19a). Hence, the fully discrete scheme (4.2) is reduced to the following:

(4.5)\begin{align} &\inf_{\boldsymbol{u}_h, \boldsymbol{u}_{T,h}} \sup_{\phi_h, \boldsymbol{\sigma}_{T,h}} \langle\!\langle\frac{|\boldsymbol{m}_h|^2}{2V_1(h_h)} +\frac{|s_h|^2}{2V_2(h_h)}\rangle\!\rangle_h + \Delta t\left\langle U(h_{T,h}) +\frac{|\boldsymbol{n}_{T,h}|^2}{2}\right\rangle_h\nonumber\\ &\qquad - \langle\!\langle h_h \partial_t\phi_h+ \boldsymbol{m}_h\boldsymbol{\cdot}\boldsymbol{\nabla}\phi_h + s_h\phi_h\rangle\!\rangle_h\nonumber\\ &\qquad -\left\langle\boldsymbol{n}_{T,h} \boldsymbol{\cdot}\boldsymbol{\sigma}_{T,h}+ \alpha h_{T,h} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma}_{T,h}\right\rangle_h\nonumber\\ &\qquad +\left\langle h_{T,h}\phi_h(T, x)-h_0\phi_h(0, x)\right\rangle_h, \end{align}

where $\boldsymbol {u}_h = (h_h, \boldsymbol {m}_h, s_h)\in [W_h^k]^{d+2}$. Moreover, the corresponding PDHG Algorithm 1 will be further simplified where the three elliptic solves in (4.4b), (4.4c) and (4.4d) are not needed, and the pointwise optimisation problem (4.4f) does not have the $\boldsymbol {n}$, $p$ and $\boldsymbol {q}$ contributions.

Since the fully discrete scheme (4.5) is first-order accurate in time, we may use a one-step approximation of the time integrals in (4.5) to reduce the computational cost without sacrificing the first-order temporal accuracy. This leads to the following approximated JKO scheme: given time step size $\Delta t >0$ and approximate solution $h_h^{n-1}\in M_h^k$ at time $t^{n-1}$, find the approximation $h_h^n\in M_h^k$ at next time level $t^n = t^{n-1}+\Delta t$ by solving the saddle-point problem

(4.6) $$\begin{align} h_h^{n}=& \arg\inf_{h_h}\inf_{\boldsymbol{u}_h} \sup_{\phi_h, \boldsymbol{\sigma}_{h}} \left\langle\frac{|\boldsymbol{m}_h|^2}{2V_1(h_h)} +\frac{|s_h|^2}{2V_2(h_h)}\right\rangle_h + \Delta t\left\langle U(h_{h}) +\frac{|\boldsymbol{n}_{h}|^2}{2}\right\rangle_h\nonumber\\ &+ \left\langle(h_h-h_h^{n-1})\phi_h - \boldsymbol{m}_h\boldsymbol{\cdot}\boldsymbol{\nabla}\phi_h - s_h\phi_h\right\rangle_h -\left\langle\boldsymbol{n}_{h} \boldsymbol{\cdot}\boldsymbol{\sigma}_{h}+ \alpha h_{h} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma}_{h}\right\rangle_h. \end{align}$$

Here $\boldsymbol {u}_h = (h_h, \boldsymbol {m}_h, s_h, \boldsymbol {n}_h)\in [M_h^k]^{2d+2}$, $\boldsymbol {\sigma }_h\in \boldsymbol {M}_{\sigma,h}^k$ and $\phi _h \in V_h^{k+1}$, where $V_h^{k+1}$ is the following $H^1$-conforming finite-element space on the spatial domain $\varOmega$: $V_h^{k+1}:=\{v\in H^1(\varOmega ): v|_{K_\ell } \in Q^{k+1}(K_\ell ) \quad \forall \ell \}.$

The saddle-point problem (4.6) can be solved using a similar optimisation solver as Algorithm 1. However, our preliminary numerical results indicate a large number of iterations (see the numerical results in § 5.1) is needed for the convergence of Algorithm 1 for solving the PDE (1.1) when the default optimisation parameters $\sigma _\phi =\sigma _u=1$ are used, especially during the time when the droplets start to merge. Here we present a direct solution strategy for (4.6) by solving the critical point system using the Newton–Raphson method. Taking the first-order variation with respect to the variables $\boldsymbol {u}_h$, $\phi _h$ and $\boldsymbol {\sigma }_h$, we get the following nonlinear system of equations: find $(h_h, \boldsymbol {m}_h, s_h, \boldsymbol {n}_h)\in [M_h^k]^{2d+2}$, $\boldsymbol {\sigma }_h\in \boldsymbol {M}_{\sigma,h}^k$ and $\phi _h\in V_h^{k+1}$ such that

(4.7a)$$\begin{gather} \left\langle\left(\frac{\boldsymbol{m}_h}{V_1(h_h)}-\boldsymbol{\nabla}\phi_h\right)\boldsymbol{\cdot}\delta\boldsymbol{m}_h\right\rangle_h=0,\quad \forall \delta\boldsymbol{m}_h\in [M_h^k]^d, \end{gather}$$

(4.7b)$$\begin{gather} \left\langle\left(\frac{s_h}{V_2(h_h)}-\phi_h\right)\boldsymbol{\cdot}\delta s_h\right\rangle_h=0,\quad \forall \delta s_h\in M_h^k, \end{gather}$$

(4.7c) $$\begin{gather}\left\langle\left(\Delta t\boldsymbol{n}_h-\boldsymbol{\sigma}_h\right)\boldsymbol{\cdot}\delta\boldsymbol{n}_h\right\rangle_h=0,\quad \forall \delta\boldsymbol{n}_h\in [M_h^k]^d, \end{gather}$$

(4.7d) $$\begin{gather} \left\langle\left(-\frac{|\boldsymbol{m}_h|^2}{2V_1^2}V_1' -\frac{|s_h|^2}{2V_2^2}V_2' +\Delta t U'(h_h)\right)\boldsymbol{\cdot}\delta h_h\right\rangle_h \nonumber\\ +\left\langle(\phi_h-\alpha\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\sigma}_h)\boldsymbol{\cdot}\delta h_h \right\rangle_h=0,\quad \forall \delta h_h\in M_h^k, \end{gather}$$

(4.7e)$$\begin{gather}\left\langle\boldsymbol{n}_{h} \boldsymbol{\cdot}\delta \boldsymbol{\sigma}_{h}+ \alpha h_{h} \boldsymbol{\nabla}\boldsymbol{\cdot}\delta\boldsymbol{\sigma}_{h}\right\rangle_h =0,\quad \forall \delta \boldsymbol{\sigma}_h\in \boldsymbol{M}_{\sigma,h}^k, \end{gather}$$

(4.7f)$$\begin{gather}\left\langle(h_h-h_h^{n-1})\delta\phi_h - \boldsymbol{m}_h\boldsymbol{\cdot}\boldsymbol{\nabla}\delta\phi_h - s_h\delta\phi_h\right\rangle_h =0,\quad \forall \delta \phi_h\in V_{h}^{k+1}. \end{gather}$$

The Newton–Raphson method is then used to solve this coupled nonlinear system.

5.1. JKO scheme for the PDE (1.1)

In this subsection, we solve the PDE (1.1) numerically using the approximated JKO scheme (4.6). Both 1-D and 2-D numerical examples are considered. In Appendix B, we present details of the spatial/temporal mesh convergence studies of the scheme (4.6).

5.1.1. 1-D example

We take the computational domain to be a periodic line segment $\varOmega =[0,1]$. The initial condition is chosen to be $h_0(x)=1-0.2\cos (2{\rm \pi} x)$, and terminal time is $T=0.4$. The approximated JKO scheme (4.6) is applied on a uniform spatial mesh with $32$ elements with polynomial degree $k=3$. In each JKO step, we use either the Newton–Raphson method to solve the critical point system (4.7) or the PDHG Algorithm 1 to solve the saddle-point problem (4.6). The PDHG iteration is terminated when the $L_1$-norm of the difference of two consecutive surface heights $h_{h}$ is less than a prescribed tolerance tol. The simulation results are compared with the following classical finite-element discretisation for (1.1) with BDF2 time stepping using the same spatial discretisation parameters: find $(h_h^n, P_h^n)\in [V_h^{k+1}]^2$ such that

(5.3a)$$\begin{gather} \left\langle\frac{3h_h^{n}-4h_h^{n-1}+h_h^{n-2}}{2\Delta t} q_h + V_1(h_h^n)\boldsymbol{\nabla} P_h^n\boldsymbol{\cdot}\boldsymbol{\nabla} q_h +V_2(h_h^n)P_h^n q_h \right\rangle_h=0, \quad \forall q_h\in V_h^{k+1}, \end{gather}$$

(5.3b)$$\begin{gather}\left\langle P_h^n r_h -U'(h_h^n)r_h-\alpha^2\boldsymbol{\nabla} h_h^n\boldsymbol{\cdot}\boldsymbol{\nabla} r_h \right\rangle_h = 0, \quad \forall r_h\in V_h^{k+1}. \end{gather}$$

In the FEM scheme (5.3), we use a small time step size $\Delta t = 0.0001$, treating these results as reference solutions. For the Newton–Raphson method, we employ both a small time step size $\Delta t = 0.0001$ and a larger time step size $\Delta t = 0.001$. In the optimisation approach, the time step size is set to $\Delta t = 0.001$, with varying stopping tolerances for each JKO iteration at $\textrm {tol} = 10^{-7}$ or $\textrm {tol} = 10^{-8}$. Figure 2 presents snapshots of surface height at different times for all these simulations. First, we observe that the results for the PDHG algorithm with $\textrm {tol}=10^{-7}$ (in magenta) differ significantly from the other simulation results, particularly at intermediate times $t=0.10$ and $t=0.15$ when the droplet is forming and spreading through dry spots. This suggests that $\textrm {tol}=10^{-7}$ is insufficient for accuracy in the PDHG algorithm. Reducing the tolerance to $\textrm {tol}=10^{-8}$ yields results consistent with the Newton–Raphson method, though the dip at the centre point $x=0.5$ at time $t=0.10$ is still not well captured compared with the reference solution. In addition, the Newton–Raphson method results with $\Delta t = 0.0001$ overlap with the FEM simulation results, validating the proposed approximate JKO scheme (4.6).

Figure 2. Snapshots of the surface height for 1-D thin-film equation (1.1) at different times: (a) $t=0.05$, (b) $t=0.10$, (c) $t=0.15$, (d) $t=0.20$, (e) $t=0.30$ and (f) $t=0.40$. Dashed black line: FEM scheme (5.3) with time step size $\Delta t = 10^{-4}$. Dotted red line: Approximated JKO scheme (4.6) with $\Delta t = 10^{-4}$ using the Newton–Raphson solver. Blue line: Approximated JKO scheme (4.6) with $\Delta t = 10^{-3}$ using the Newton–Raphson solver. Green line: Approximated JKO scheme (4.6) with $\Delta t = 10^{-3}$ using the PDHG solver with tolerance $\textrm {tol}=10^{-8}$. Magenta line: Approximated JKO scheme (4.6) with $\Delta t = 10^{-3}$ using the PDHG solver with tolerance $\textrm {tol}=10^{-7}$.

5.1.2. 2-D example

We consider the computational domain as a periodic unit square $\varOmega =[0,1]^2$. The initial condition is chosen as $h_0(x,y)=1+0.2\cos (2{\rm \pi} x)\cos (2{\rm \pi} y)$, and the terminal time is $T=0.4$. The approximated JKO scheme (4.6) is applied on a uniform rectangular mesh with $32\times 32$ elements and polynomial degree $k=3$. Similar to the 1-D case, we use either the Newton–Raphson method to solve the critical point system (4.7) or the PDHG Algorithm 1 to solve the saddle-point problem (4.6) for each JKO step.

We compare the simulation results with the reference solution using the FEM (5.3) with the same spatial discretisation parameters. For the FEM scheme (5.3), we use a small time step size $\Delta t = 0.0001$, and these results are treated as reference solutions. For the Newton–Raphson method, we use both a small time step size $\Delta t = 0.0001$ and a large time step size $\Delta t = 0.001$. For the optimisation approach, we set the time step size to be $\Delta t = 0.001$, and terminate the iteration when the $L_1$-norm of the difference of two consecutive surface heights $h_{h}$ is less than a tolerance $\textrm {tol}=10^{-8}$.

Snapshots of surface height contours at different times are shown in figure 3. We observe a similar pattern for the surface height evolution as in the 1-D case. Driven by the interfacial instabilities described in (1.1), the early stage (see figure 3a) shows the spatial variations in the solution profile growing until the minimum height approaches $h_{min} = O(\epsilon )$. In the later stage, the minimum height spreads to form dry spots, leading to droplet formation (see figure 3b), followed by a slow growth in droplet height driven by weak condensation effects (see figure 3c–f). This example captures the morphological changes previously observed in 1-D dewetting thin-film dynamics with weak non-mass-conserving effects (Ji & Witelski Reference Ji and Witelski2018). Moreover, results for all four simulations are qualitatively similar to each other. In particular, the results for the first two rows with the small time step size $\Delta t=0.0001$ are almost identical to each other, validating the proposed approximated JKO scheme (4.6) in the 2-D setting. The results for the last two rows are consistent with each other, which indicates both the Newton–Raphson and the PDHG optimisation approaches are able to solve the discrete saddle-point problem (4.6) accurately.

Figure 3. Snapshots of the surface height contour for the 2-D thin-film equation (1.1) at different times: (a) $t=0.05$, (b) $t=0.10$, (c) $t=0.15$, (d) $t=0.20$, (e) $t=0.30$ and ( f) $t=0.4$; (i) numerical solutions for FEM (5.3), $\Delta t=0.0001$; (ii) numerical solutions for the approximated JKO scheme (4.6) with the Newton–Raphson solver for (4.7) in each JKO step, $\Delta t=0.0001$; (iii) numerical solutions for the approximated JKO scheme (4.6) with the Newton–Raphson solver for (4.7) in each JKO step, $\Delta t=0.001$; (iv) numerical solutions for the approximated JKO scheme (4.6) with PDHG solver in Algorithm 1 for each JKO step, $\Delta t=0.001$, $\textrm {tol} = 10^{-8}$.

5.2. MFC for droplet dynamics

Next, we demonstrate that the developed MFC system can drastically control droplet shapes, motions and drive morphological changes in droplet configurations. With different choices of initial and target surface heights, we present numerical results showcasing the application of the developed MFC scheme for fundamental droplet actuation techniques, including droplet transport, bead-up (i.e. dewetting), spreading, merging and splitting. To achieve this, we apply the finite-element discretisation (4.2) and the PDHG Algorithm 1 to solve the full MFC system in Definition 3.3. We use the same mobility functions and energy functional (5.1) as in the previous subsection. There is additional freedom in choosing the potential and terminal functionals in (3.4a) and the scaling constant $\beta$ in (3.4b) to determine the MFC problem (3.4). Throughout this subsection, we take the terminal time $T=1$ and $\beta =0.01$, and use the following potential and terminal functionals:

(5.4)\begin{equation} \mathcal{F}(h) ={-}0.02\int_0^T\int_{\varOmega}h\log(h) \,{{\rm d}\kern0.7pt x}\,{\rm d}t,\quad \mathcal{G}(h) = 0.05\int_{\varOmega}h(\log(h/{h_{T}})-1) \,{{\rm d}\kern0.7pt x}, \end{equation}

where ${h_T}$ is a given target terminal surface height function. The potential functional $\mathcal {F}(h)$ serves as a regularisation term, whereas the terminal function $\mathcal {G}(h)$ drives the surface height towards the target surface height as $h_T=\arg \min _h \mathcal {G}(h)$. The choice of these functional forms for different control constraints is beyond the scope of this study and will serve as an interesting future direction.

In the finite-element discretisation (4.2), we use a uniform spatial rectangular mesh of size $64\times 64$, a uniform temporal mesh of size $16$, and take polynomial degree $k=3$. The total number of degrees of freedom (DOFs) for the dual variable $\phi _h$ is about 4 million. This choice of discretisation parameters ensures enough resolution in the simulation with a manageable computational cost. A more in-depth mesh resolution study will be carried out elsewhere. We stop the PDHG iteration when the $L_1$-norm of the difference between the terminal surface heights in two consecutive iterations $\|h_{h,T}^n-h_{h,T}^{n-1}\|_{L_1}$ is less than $10^{-6}$, when a converged solution profile has been reached. For all the test cases considered here, the algorithm converged within 2000 iterations. As mentioned in the previous subsection, a more in-depth computational study of the optimisation solver Algorithm 1 will be carried elsewhere.

5.2.1. Case 1: droplet transport

Droplet transport is one of the most important operations in DMF and has been investigated extensively through experimental approaches such as electro-dewetting (Li et al. Reference Li, Ha, Liu, van Dam and Kim2019). We present results for the MFC problem (3.4) with the initial and target surface heights specified as

(5.5)$$\begin{gather} h_0(x, y) = \epsilon+ \frac{10}{3}\left(1-75((x-0.3)^2+(y-0.3)^2)\right)_+, \end{gather}$$

(5.6)$$\begin{gather}{h_{T}(x, y)} = \epsilon+ \frac{10}{3}\left(1-75((x-0.7)^2+(y-0.7)^2)\right)_+. \end{gather}$$

Here

(5.7) \begin{equation} (f)_+= \begin{cases} f, & \text{if } f\ge0,\\ 0, & \text{otherwise}, \end{cases} \end{equation}

is the positive part of a function $f$. This example models the MFC of an initial parabolic droplet centred at $(0.3, 0.3)$ moving towards a target parabolic droplet centred at $(0.7,0.7)$. Figure 4 shows the evolution of the controlled surface height at different times, alongside the activity field $\zeta _h=({\beta (U'(h_h)+p_h)+\phi _h})/{h_h}$ in (3.23). The snapshots in figure 4 reveal that the droplet quickly breaks the symmetry, beads up and develops a larger advancing contact angle, leaving a capillary wave as it progresses towards the target position. The three-dimensional (3-D) contour plots in figure 4(i) are reconstructed using piecewise polynomial interpolation of the discrete surface height data in the space–time quadrature space $h_h\in W_h^k$. Specifically, a third-order polynomial interpolation ($k=3$) is employed for the surface reconstruction within each space–time cell, which is then used to evaluate the surface height at specific times. Due to the non-smoothness of the initial profile (5.5), visible oscillations in the reconstructed surface height are observed near the contact line (where strong gradients exist) at the initial time $t=0$. These oscillations are localised near the contact lines and diminish as time progresses due to the instantaneous smoothing effect of the nonlinear diffusion terms in the MFC problem. In addition, the 2-D contour plots in figures 4(ii) and 4(iii) are reconstructed from piecewise constant interpolation of the data on a refined $256\times 256$ mesh. Similar to the 3-D plots, localised oscillations near the contact lines are observed, which diminish over time. The surface height remains positive throughout the simulation, achieving a smooth profile by the final time $t=1$. Furthermore, the activity field tends to be contractile ($\zeta > 0$) near the target droplet profile and extensile ($\zeta < 0$) near the initial profile, as the active stress competes with passive interfacial forces and pushes the droplet towards the terminal location.

Figure 4. Case 1 (droplet transport) snapshots of (i) 3-D plots of the controlled surface height $h$, (ii) contour plots of the controlled surface height and (iii) contour plots of the activity field $\zeta$ at different times: (a) $t=0.0$, (b) $t=0.2$, (c) $t=0.4$, (d) $t=0.6$, (e) $t=0.8$ and (f) $t=1.0$. The corresponding animation video can be found in the GitHub repository (Fu et al. Reference Fu, Ji, Pazner and Li2024a)

5.2.2. Case 2: droplet spreading

Controlling the deformation of a droplet is also a crucial aspect of liquid-handling technology. For instance, in typical electro-wetting and electro-dewetting experiments, an electric field can induce changes in the contact angles of a slender droplet containing a dilute surfactant (Nelson & Kim Reference Nelson and Kim2012). Here, we demonstrate the MFC mechanism to control the spreading and bead-up of droplets.

For the droplet-spreading example, we take the initial and target surface heights as

(5.8)$$\begin{gather} h_0(x, y) = \epsilon+\frac{10}{3}\left(1-75((x-0.5)^2+(y-0.5)^2)\right)_+, \end{gather}$$

(5.9)$$\begin{gather}{h_{T}(x, y)} = \epsilon+ \frac{5}{12}\left(1-\frac{75}{8}((x-0.5)^2+(y-0.5)^2)\right)_+. \end{gather}$$

This case models the MFC of an initial parabolic droplet centred at $(0.5, 0.5)$ with half-width $w=\tfrac {1}{5\sqrt {3}}$ flattening towards the target droplet with half-width $w=\tfrac {2\sqrt {2}}{5\sqrt {3}}$. The initial and target droplets have the same total mass. Snapshots of the simulation results for the scheme (4.2) are presented in figure 5. The numerical results indicate that the controlled droplet initially evolves into a pancake shape and then gradually converges to the target droplet profile over time. The radial symmetry in the droplet profile is preserved during the evolution. Compared with the previous droplet transport test case, visible oscillations are only observed at initial time $t=0$, as shown in figure 4(a i). In this test case, the activity field remains radially symmetric and contractile ($\zeta > 0$) throughout the simulation, reaching its maximum value near the droplet's contact lines at each time.

Figure 5. Case 2 (droplet spreading) snapshots of (i) 3-D plots of the controlled surface height $h$, (ii) contour plots of the controlled surface height and (iii) contour plots of the activity field $\zeta$ at different times: (a) $t=0.0$, (b) $t=0.2$, (c) $t=0.4$, (d) $t=0.6$, (e) $t=0.8$ and (f) $t=1.0$. The corresponding animation video can be found in the GitHub repository (Fu et al. Reference Fu, Ji, Pazner and Li2024a).

5.2.3. Case 3: droplet bead-up

This represents the reverse process of case 2, where our objective is to induce droplet bead-up (i.e. dewetting). Therefore, we set the initial and target surface heights as

(5.10)$$\begin{gather} h_0(x, y) = \epsilon+ \frac{5}{12}\left(1-\frac{75}{8}((x-0.5)^2+(y-0.5)^2)\right)_+, \end{gather}$$

(5.11)$$\begin{gather}{h_{T}(x, y)} = \epsilon+ \frac{10}{3}\left(1-75((x-0.5)^2+(y-0.5)^2)\right)_+, \end{gather}$$

where the initial and target profiles are reversed compared with the example in case 2. This models the MFC of an initial parabolic droplet centred at $(0.5, 0.5)$ with half-width $w= \tfrac {2\sqrt {2}}{5\sqrt {3}}$ which beads up and evolves into the target droplet with half-width $w=\tfrac {1}{5\sqrt {3}}$. Snapshots of the simulation results for the scheme (4.2) are presented in figure 6. In this case, we observe more pattern formation during the droplet bead-up process, where capillary waves are generated in the early stage before the droplet reaches the target shape. Compared with the target droplet profile, the obtained profile at the terminal time $t = 1$ has a slightly elevated height near the contact line. Moreover, the activity field remains radially symmetric and mostly extensile ($\zeta < 0$) throughout the simulation, reaching its maximum absolute value near the droplet's contact lines at each time.

Figure 6. Case 3 (droplet bead-up) snapshots of (i) 3-D plots of the controlled surface height $h$, (ii) contour plots of the controlled surface height and (iii) contour plots of the activity field $\zeta$ at different times: (a) $t=0.0$, (b) $t=0.2$, (c) $t=0.4$, (d) $t=0.6$, (e) $t=0.8$ and (f) $t=1.0$. The corresponding animation video can be found in the GitHub repository (Fu et al. Reference Fu, Ji, Pazner and Li2024a).

5.2.4. Case 4: droplet merging

Finally, we demonstrate the application of MFC for droplet merging and splitting, which are more complex droplet manipulation techniques widely used in biological and chemical applications (Nan, Mao & Shum Reference Nan, Mao and Shum2023). For the droplet-merging case, we control the coalescence of four small droplets initially placed on a 2-D domain, where the initial surface height profile is

(5.12a)\begin{equation} h_0(x, y) = \epsilon+\sum_{i=1}^{4}\tfrac{5}{6}[1-75((x-x_i)^2+(y-y_i)^2)]_+, \end{equation}

where the positions of the peaks of the initial droplets are $(x_1,y_1) = (0.3, 0.3)$, $(x_2,y_2) = (0.3, 0.7)$, $(x_3, y_3) = (0.7,0.3)$ and $(x_4, y_4) = (0.7,0.7)$.

We consider two subcases for the terminal surface height: the first subcase has a symmetric terminal droplet profile centred at $(0.5, 0.5)$ which has the same total mass as the initial profile,

(5.12b)\begin{equation} h_{T}(x, y) = \epsilon+ \frac{5}{12}\left(1-\frac{75}{8}((x-0.5)^2+(y-0.5)^2)\right)_+, \end{equation}

whereas the second subcase has a skewed droplet profile elongated in the $x$-direction centred at $(0.6,0.55)$, which has a different total mass as the initial profile,

(5.12c)\begin{equation} h_{T}(x, y) = \epsilon+ \frac{5}{12}\left(1-\frac{75}{8}(4(x-0.6)^2+(y-0.55)^2)\right)_+, \end{equation}

The second subcase showcases the effect of asymmetry in the evolution of the controlled surface height.

Figure 7 presents snapshots of the simulation results for the first subcase with the symmetric target surface height (5.12b). Similar to the droplet transport example discussed in case 1, the controlled droplets quickly bead up and start moving towards the centre of the domain, where the target droplet is placed. Capillary waves form behind the droplets as they shift towards the target position. Near the terminal time $t = 1$, the droplets coalesce, forming a single droplet at the centre of the domain. Notably, in this test case, contact line oscillations are only observed at the initial time $t=0$ due to the non-smoothness of the initial data. A radially symmetric and smooth solution profile for $h(t, \boldsymbol {x})$ is observed for $t > 0$. The activity field $\zeta$ at different times is plotted in figure 7(iii). This shows that in the early stage, high contractile active stress ($\zeta > 0$) appears near the centre of the domain, whereas extensile active stress ($\zeta < 0$) around the individual droplets pushes the fluid towards the centre of the domain. As the droplets become closer to each other, the corresponding contractile active stress becomes more concentrated near the centre of the domain, and the magnitude of the activity field diminishes as the terminal droplet is formed.

Figure 7. Case 4 (droplet merging with a symmetric target profile (5.12b)) snapshots of (i) 3-D plots of the controlled surface height $h$, (ii) contour plots of the controlled surface height and (iii) contour plots of the activity field $\zeta$ at different times: (a) $t=0.0$, (b) $t=0.2$, (c) $t=0.4$, (d) $t=0.6$, (e) $t=0.8$ and (f) $t=1.0$. The corresponding animation video can be found in the GitHub repository (Fu et al. Reference Fu, Ji, Pazner and Li2024a).

Figure 8 presents simulation snapshots for the second subcase with the asymmetric target surface height (5.12c). This case exhibits more interesting transient dynamics due to the asymmetry in the target droplet located at the off-centre position $(x, y) = (0.6, 0.55)$. In the early stage, the two droplets in the left half of the domain quickly shrink in size, whereas the two droplets in the right half of the domain gradually bead up at different speeds. In the later stage, the two remaining droplets slowly move towards each other and coalesce into the target droplet, which is elongated in the $x$ direction at the desired off-centre location at time $t = 1$. The corresponding activity field plots illustrate that the high extensile active stress ($\zeta < 0$) around the two droplets on the left contributes to the rapid vanishing of the droplets. Meanwhile, the combined effects of the high contractile active stress ($\zeta > 0$) near the location of the terminal droplet and the low extensile stress ($\zeta < 0$) near the contact line of the two droplets on the right lead to the directed motion of the droplets towards the terminal location. This example specifically showcases the flexibility of our control mechanism in handling initial and terminal configurations of different total masses through the non-mass-conserving effects.

Figure 8. Case 4 (droplet merging with an asymmetric target profile (5.12c)) snapshots of (i) 3-D plots of the controlled surface height $h$, (ii) contour plots of the controlled surface height and (iii) contour plots of the activity field $\zeta$ at different times: (a) $t=0.0$, (b) $t=0.2$, (c) $t=0.4$, (d) $t=0.6$, (e) $t=0.8$ and (f) $t=1.0$. The corresponding animation video can be found in the GitHub repository (Fu et al. Reference Fu, Ji, Pazner and Li2024a).

5.2.5. Case 5: droplet splitting

For the case of droplet splitting, we consider the reverse process of the symmetric droplet merging considered in case 4 and take the initial and target surface heights as

(5.13)$$\begin{gather} h_0(x, y) = \epsilon+ \frac{5}{12}\left[1-\frac{75}{8}((x-0.5)^2+(y-0.5)^2)\right]_+, \end{gather}$$

(5.14)$$\begin{gather}{h_{T}(x, y)} = \epsilon+\sum_{i=1}^{4}\frac{5}{6}[1-75((x-x_i)^2+(y-y_i)^2)]_+, \end{gather}$$

where the initial and target profiles are reversed compared with the first example in case 4, with the droplet configurations in the initial and target profiles identical to those specified in (5.12c) and (5.12a), respectively. Snapshots of the simulation results for the controlled dynamics are presented in figure 9. This example models the MFC of one initial parabolic droplet splitting into four smaller droplets and shifting into target positions individually. Compared with the droplet-merging case, splitting a single droplet appears to be more challenging and the terminal surface height profile obtained at the terminal time $t = 1$ still maintains a ridge connecting the small droplets. The activity plots in figure 9(iii) demonstrate that near time $t = 0$, the inhomogeneous distribution of high extensile active stress ($\zeta < 0$) near the contact line of the centre droplet contributes to the fast splitting of the initial droplet. Meanwhile, the high contractile active stress ($\zeta > 0$) around the terminal locations of the target droplets drives the individual droplets away from the centre and pushes them towards the target locations. As the terminal time $t = 1$ approaches, high extensile stress ($\zeta < 0$) at the centre of the domain and high contractile stress ($\zeta > 0$) are formed as the individual droplets spread and settle into their target shapes.

Figure 9. Case 5 (droplet splitting) snapshots of (i) 3-D plots of the controlled surface height $h$, (ii) contour plots of the controlled surface height and (iii) contour plots of the activity field $\zeta$ at different times: (a) $t=0.0$, (b) $t=0.2$, (c) $t=0.4$, (d) $t=0.6$, (e) $t=0.8$ and (f) $t=1.0$. The corresponding animation video can be found in the GitHub repository (Fu et al. Reference Fu, Ji, Pazner and Li2024a); see supplementary movie available at https://doi.org/10.1017/jfm.2024.983.