Phase-field and dissipative gradient-flow computation

Let $\mathrm{\Omega }\subset {\mathbb{R}}^d$ be a bounded connected Lipschitz domain, let $d\ge 1$, and fix constants $\epsilon >0$, $M>0$, and $\kappa >0$. Put $F_{\epsilon }(s)=\frac{(s^2-1{)}^2}{4{\epsilon }^2},f_{\epsilon }(s)=F_{\epsilon }^{\mathrm{\prime }}(s)=\frac{s^3-s}{{\epsilon }^2},$ and define the free energy on $H^1(\mathrm{\Omega };\mathbb{R})\cap L^4(\mathrm{\Omega };\mathbb{R})$ by ${\mathcal{E}}_{\epsilon }(u)={\int }_{\mathrm{\Omega }}(\frac{\kappa }{2}\mathrm{\mid }\mathrm{\nabla }u{\mathrm{\mid }}^2+F_{\epsilon }(u))dx\mathrm{.}$ Here $H^1=W^{1,2}$ and weak derivatives have the meaning in Integer-order Sobolev spaces and their norms. The finite-energy initial class consists of $u_0\in H^1(\mathrm{\Omega })\cap L^4(\mathrm{\Omega })$.

The Allen–Cahn gradient flow is $u_t=-M\mu ,\mu =-\kappa \mathrm{\Delta }u+f_{\epsilon }(u),{\mathrm{\partial }}_{n}u=0,u(0)=u_0\mathrm{.}$ The Cahn–Hilliard gradient flow is $u_t=M\mathrm{\Delta }\mu ,\mu =-\kappa \mathrm{\Delta }u+f_{\epsilon }(u),{\mathrm{\partial }}_{n}u={\mathrm{\partial }}_n\mu =0,u(0)=u_0\mathrm{.}$ For weak solutions, the no-flux conditions are encoded by the natural weak formulations with all $H^1(\mathrm{\Omega })$ test functions. In dimensions $d\ge 5$, the finite-energy class $H^1\cap L^4$ alone does not in general make $f_{\epsilon }(u)$ a continuous functional on $H^1(\mathrm{\Omega })$: it gives only $f_{\epsilon }(u)\in L^{4\mathrm{/}3}$, while $H^1(\mathrm{\Omega })$ need not embed into $L^4(\mathrm{\Omega })$. Thus any weak solution tested against all $H^1$ functions must have the additional integrability needed to make the nonlinear pairing continuous. The chemical potential $\mu$ is the variational derivative of ${\mathcal{E}}_{\epsilon }$, and the mobilities select the $L^2$ metric for Allen–Cahn and the mean-zero $H^{-1}$ metric for Cahn–Hilliard. In the Cahn–Hilliard flow, the no-flux equation conserves ${\int }_{\mathrm{\Omega }}u$; the Allen–Cahn flow has no such conservation law in general. No maximum-principle property is part of these definitions.

The convergence theorem on this page is restricted to one space dimension. Assuming the Axiom of Choice for the following standard embedding consequences, if $d=2$ and $\mathrm{\Omega }$ is a bounded $W^{1,p}$ extension domain for every $1<p<2$ as defined in Sobolev extension domains and extension operators, then The critical Sobolev embedding into every finite $L^q$ gives $H^1(\mathrm{\Omega })\hookrightarrow L^4(\mathrm{\Omega })$. If $d=3$ and $\mathrm{\Omega }$ is a bounded $W^{1,2}$ extension domain, then Sobolev embedding on bounded extension domains for $p<n$ gives $H^1(\mathrm{\Omega })\hookrightarrow L^4(\mathrm{\Omega })$. These are the exact extra extension hypotheses for the quartic-energy claim. In dimensions $d\ge 4$ no assertion that every $H^1$ state has finite quartic energy is made; the explicit initial class above remains $H^1\cap L^4$. No sharp-interface limit as $\epsilon \downarrow 0$ is part of either definition.

The no-flux finite-volume realization starts from a finite connected cell-face mesh as in Finite-volume cell-face meshes and geometric regularity and uses the interior-face orthogonality and positive transmissibilities of Orthogonal TPFA diffusion with harmonic face transmissibility. For unit diffusivity, let $X_h={\mathbb{R}}^{\mathcal{K}}$, $m_h(u,v)={\sum }_{K\in \mathcal{K}}\mathrm{\mid }K\mathrm{\mid }{u}_K{v}_K$, and $a_h(u,v)={\sum }_{F=K\mathrm{\mid }L\in {\mathcal{F}}_{\mathrm{i}\mathrm{n}\mathrm{t}}}{T}_F(u_K-u_L)(v_K-v_L)\mathrm{.}$ There is no boundary-face term: the flux through each physical boundary face is prescribed to be zero. Define $A_h:X_h\to X_h$ by $m_h(A_{h}u,v)=a_h(u,v)$ for all $v\in X_h$. Split the double-well density as $F_{\epsilon }=F_c-F_e,F_c(s)=\frac{s^4}{4{\epsilon }^2}+\frac{s^2}{2{\epsilon }^2},F_e(s)=\frac{s^2}{{\epsilon }^2}-\frac{1}{4{\epsilon }^2}\mathrm{.}$ The associated finite-volume energy components are $E_{c,h}(u)={\sum }_K\mathrm{\mid }K\mathrm{\mid }(\frac{u_K^4}{4{\epsilon }^2}+\frac{u_K^2}{2{\epsilon }^2}),E_{e,h}(u)={\sum }_K\mathrm{\mid }K\mathrm{\mid }(\frac{u_K^2}{{\epsilon }^2}-\frac{1}{4{\epsilon }^2})\mathrm{.}$

For convergence on $\mathrm{\Omega }=(0,L)$, a second realization uses the same orthogonal dual mesh with centers $x_i=(i-\frac{1}{2})h$, $h=L\mathrm{/}N$. Its conforming space $V_h$ consists of functions affine between consecutive centers and constant on the two half-cells between a boundary and the nearest center. Its inner product is the exact $L^2(0,L)$ inner product, its stiffness is $a_h(v,w)={\int }_0^L{v}^{\mathrm{\prime }}{w}^{\mathrm{\prime }}dx$, and its energy components are $E_{c,h}(v)={\int }_0^L{F}_c(v)dx$ and $E_{e,h}(v)={\int }_0^L{F}_e(v)dx$. Define $A_h:V_h\to V_h$ by $(A_{h}v,w{)}_{L^2}=a_h(v,w)$. For degrees of freedom $v_i=v(x_i)$, this stiffness is exactly the interior two-point flux form ${\sum }_{i=1}^{N-1}(v_{i+1}-v_i)(w_{i+1}-w_i)\mathrm{/}h$; the exact mass and potential integrals are the conforming reconstruction used in the convergence theorem.

Both realizations have a constant function $\mathbf{1}$, a real finite-dimensional inner-product space, and a self-adjoint nonnegative operator $A_h$. They share the physical energy $E_h(u)=\frac{\kappa }{2}{m}_h(A_{h}u,u)+E_{c,h}(u)-E_{e,h}(u)$, with the appropriate inner product and energy components for the realization. The sign convention is $A_h\approx -\mathrm{\Delta }$.

The diagram labels four adjacent cell values, three shared interior fluxes, and the prescribed zero boundary fluxes. Each interior flux appears with opposite signs in its two cell balances. The diagram records the convention only; it is not used as a proof.

The published TPFA theorem TPFA discrete Green identity, coercivity, and maximum principle treats homogeneous Dirichlet data and its energy includes boundary penalties. The no-flux realization above deliberately omits those terms and uses only its interior transmissibilities; its kernel and gauge are proved separately in Kernel, range, and mean-zero inverse for the no-flux Laplacian. The conditional parabolic limit item Conditional monotone finite-volume parabolic limit identification assumes zero-boundary $W_0^{1,p}$ reconstructions and is not used for this Neumann branch. The finite-dimensional gradient-flow algorithm also applies to the conforming realization stated above.

Use either connected no-flux realization in No-flux finite-volume and conforming realizations for phase fields. In its finite-volume branch, assume the interior cell graph is connected; in its interval conforming branch, assume $N\ge 1$. Write ${\bar{u}}_h=⟨u_h,\mathbf{1}{⟩}_h\mathrm{/}⟨\mathbf{1},\mathbf{1}{⟩}_h$ and $X_{h,0}=\{u_h:{\bar{u}}_h=0\}$. Then $A_h$ is self-adjoint and nonnegative, $\ker A_h=\mathrm{span}\{\mathbf{1}\}$, and $\mathrm{ran}{A}_h=X_{h,0}$. The restriction $A_h:X_{h,0}\to X_{h,0}$ is positive definite and invertible. Define $A_h^{†}{u}_h$ by first projecting $u_h$ onto $X_{h,0}$ and then applying this self-adjoint positive inverse. For each fixed finite mesh there is $C_h<\mathrm{\infty }$ such that $\mathrm{\parallel }{u}_h-{\bar{u}}_h\mathbf{1}{\mathrm{\parallel }}_h^2\le C_h⟨A_h{u}_h,u_h{⟩}_h\mathrm{.}$ For a one-cell mesh, $X_{h,0}=\{0\}$ and the inequality is vacuous.

Let $\mathrm{\Omega }=(0,L)$ and let $V_h$ be the conforming interval space in No-flux finite-volume and conforming realizations for phase fields, with center spacing $h$. Let $I_{h}v$ interpolate the absolutely-continuous representative of $v\in H^1(0,L)$ at the centers and extend it constantly to the two boundary half-cells. Let $P_h$ be the exact $L^2$-orthogonal projection onto $V_h$. There is a constant $C$ independent of $h$ such that $\mathrm{\parallel }{I}_{h}v{\mathrm{\parallel }}_{H^1}\le C\mathrm{\parallel }v{\mathrm{\parallel }}_{H^1},\mathrm{\parallel }v-I_{h}v{\mathrm{\parallel }}_{L^2}\le Ch\mathrm{\parallel }{v}^{\mathrm{\prime }}{\mathrm{\parallel }}_{L^2},\mathrm{\parallel }{P}_{h}v{\mathrm{\parallel }}_{H^1}\le C\mathrm{\parallel }v{\mathrm{\parallel }}_{H^1}\mathrm{.}$ In addition, $\mathrm{\parallel }{P}_{h}v-v{\mathrm{\parallel }}_{L^2}\le Ch\mathrm{\parallel }{v}^{\mathrm{\prime }}{\mathrm{\parallel }}_{L^2}$ for every $v\in H^1(0,L)$, $P_{h}v\to v$ in $H^1(0,L)$ for every $v\in H^1(0,L)$, and $P_{h}g\to g$ in $L^2(0,L)$ for every $g\in L^2(0,L)$ as $h\to 0$.

On either discrete realization of No-flux finite-volume and conforming realizations for phase fields, write $X_h={\mathbb{R}}^{\mathcal{K}}$ for the finite-volume state space and $X_h=V_h$ for the conforming interval state space; use its inner product $⟨\cdot ,\cdot {⟩}_h$ and the energy components $E_{c,h},E_{e,h}$ defined there. Let $g_h(u)={\mathrm{\nabla }}_h(E_{c,h}-E_{e,h})(u)$ be the Riesz gradient of the full potential energy. In the finite-volume realization, $(g_h(u){)}_K=f_{\epsilon }(u_K)$. In the conforming interval realization, $g_h(u)\in V_h$ is defined by $⟨g_h(u),v_h{⟩}_h={\int }_0^L{f}_{\epsilon }(u_h(x))v_h(x)dx(v_h\in V_h),$ where $⟨\cdot ,\cdot {⟩}_h$ is the exact $L^2$ inner product. Given a previous state $u^n\in X_h$, parameters $\epsilon ,\kappa ,M,\tau >0$, and a stabilization parameter $S\ge 0$, the stabilized semi-implicit Allen–Cahn step is the linear equation $\frac{u^{n+1}-u^n}{\tau }=-M(\kappa A_h{u}^{n+1}+g_h(u^n)+S(u^{n+1}-u^n))\mathrm{.}$ For the same input, the stabilized semi-implicit Cahn–Hilliard step is $\frac{u^{n+1}-u^n}{\tau }=-M{A}_h{\mu }^{n+1},{\mu }^{n+1}=\kappa A_h{u}^{n+1}+g_h(u^n)+S(u^{n+1}-u^n)\mathrm{.}$ These formulas define alternative updates with an explicitly lagged nonlinear potential and an added increment term. No energy law, unconditional solvability, maximum principle, or error estimate is asserted here for the quartic double well; any such property requires its own hypotheses and proof. In particular, no stability estimate proved for the convex-splitting update is silently transferred to this variant.

Use either Cahn–Hilliard realization of No-flux finite-volume and conforming realizations for phase fields. Write $X_h={\mathbb{R}}^{\mathcal{K}}$ for the finite-volume state space and $X_h=V_h$ for the conforming interval state space; use the specified inner product $⟨\cdot ,\cdot {⟩}_h$ in each case. Set $E_{1,h}=E_{c,h}-E_{e,h}$. Since $F_{\epsilon }(s)=(s^2-1{)}^2\mathrm{/}(4{\epsilon }^2)\ge 0$, this is nonnegative. Explicitly, $E_{1,h}(u)={\sum }_{K\in \mathcal{K}}\mathrm{\mid }K\mathrm{\mid }{F}_{\epsilon }(u_K)\text{in the finite-volume realization},$ and $E_{1,h}(u_h)={\int }_0^L{F}_{\epsilon }(u_h(x))dx\text{in the conforming interval realization}\mathrm{.}$ Its Riesz gradient $g_{1,h}(u)={\mathrm{\nabla }}_h{E}_{1,h}(u)$ is specified in the finite-volume realization by $(g_{1,h}(u){)}_K=f_{\epsilon }(u_K)$, where $f_{\epsilon }=F_{\epsilon }^{\mathrm{\prime }}$. In the conforming interval realization it is the unique $g_{1,h}(u_h)\in V_h$ satisfying $⟨g_{1,h}(u),v_h{⟩}_h={\int }_0^L{f}_{\epsilon }(u_h(x))v_h(x)dx(v_h\in V_h),$ with the exact $L^2$ inner product from that realization.

Inputs: $u^n\in X_h$, $r^n\in \mathbb{R}$, and parameters $\epsilon ,\kappa ,M,\tau ,C_0>0$. Set ${\alpha }_n=(E_{1,h}(u^n)+C_0{)}^{-1\mathrm{/}2}$. Output: the tuple $(u^{n+1},r^{n+1},{\mu }^{n+1})$ obtained by solving the linear equations $\frac{u^{n+1}-u^n}{\tau }=-M{A}_h{\mu }^{n+1},$ ${\mu }^{n+1}=\kappa A_h{u}^{n+1}+{\alpha }_n{r}^{n+1}{g}_{1,h}(u^n),$ $r^{n+1}-r^n=\frac{{\alpha }_n}{2}⟨g_{1,h}(u^n),u^{n+1}-u^n{⟩}_h\mathrm{.}$ The initial scalar may be set to $r^0=\sqrt{E_{1,h}(u^0)+C_0}$. After a step, $r^{n+1}$ is an independent auxiliary variable determined by the third equation; the scheme does not impose $r^{n+1}=\sqrt{E_{1,h}(u^{n+1})+C_0}$.

The separate result The SAV step is uniquely solvable and dissipates a modified energy ↗ proves unique solvability of this linear system and its modified-energy identity.

Let $X_h$ be either finite-dimensional realization in No-flux finite-volume and conforming realizations for phase fields, with real inner product $⟨\cdot ,\cdot {⟩}_h$, constant vector $\mathbf{1}$, self-adjoint nonnegative $A_h$ with kernel $\mathrm{span}\{\mathbf{1}\}$, and $A_h^{†}$ from Kernel, range, and mean-zero inverse for the no-flux Laplacian. Let $P_{c,h},P_{e,h}:X_h\to \mathbb{R}$ be the differentiable convex energy components specified there; ${\mathrm{\nabla }}_{h}P$ denotes the Riesz gradient relative to $⟨\cdot ,\cdot {⟩}_h$. Put ${\mathcal{E}}_h(u)=\frac{\kappa }{2}⟨A_{h}u,u{⟩}_h+P_{c,h}(u)-P_{e,h}(u),{\mu }_h(u,v)=\kappa A_{h}u+{\mathrm{\nabla }}_h{P}_{c,h}(u)-{\mathrm{\nabla }}_h{P}_{e,h}(v)\mathrm{.}$ Given $u^n\in X_h$, $\tau >0$, and $M>0$, the backward-Euler convex-splitting step is:

1. Allen–Cahn branch. Return the minimizer over $v\in X_h$ of $J_{AC}(v)=\frac{1}{2M\tau }\mathrm{\parallel }v-u^n{\mathrm{\parallel }}_h^2+\frac{\kappa }{2}⟨A_{h}v,v{⟩}_h+P_{c,h}(v)-⟨{\mathrm{\nabla }}_h{P}_{e,h}(u^n),v{⟩}_h,$ and set ${\mu }^{n+1}={\mu }_h(u^{n+1},u^n)$.

2. Cahn–Hilliard branch. Restrict to the affine set $⟨v,\mathbf{1}{⟩}_h=⟨u^n,\mathbf{1}{⟩}_h$ and return the minimizer of $J_{CH}(v)=\frac{1}{2M\tau }⟨A_h^{†}(v-u^n),v-u^n{⟩}_h+\frac{\kappa }{2}⟨A_{h}v,v{⟩}_h+P_{c,h}(v)-⟨{\mathrm{\nabla }}_h{P}_{e,h}(u^n),v{⟩}_h,$ then set ${\mu }^{n+1}={\mu }_h(u^{n+1},u^n)$.

For the finite-volume realization, $P_{c,h}(u)={\sum }_K\mathrm{\mid }K\mathrm{\mid }{F}_c(u_K)$ and $P_{e,h}(u)={\sum }_K\mathrm{\mid }K\mathrm{\mid }{F}_e(u_K)$. For the conforming interval realization, these are the exact integrals of $F_c$ and $F_e$ over the reconstructed piecewise-affine state. The recurrence, minimizer equivalence, uniqueness, conservation, and energy bounds are proved in Existence, uniqueness, mass conservation, and energy decay of a convex-splitting step ↗. This is the selected backward-Euler branch; it asserts no ARS(2,2,2), Lie, or Strang property.

Under the finite-dimensional hypotheses of First-order scalar-auxiliary-variable Cahn–Hilliard step, the displayed SAV linear system has a unique solution at every step. Every solution conserves the discrete mean and obeys ${\tilde{E}}_h(u^{n+1},r^{n+1})+\frac{\kappa }{2}⟨A_h\delta u,\delta u{⟩}_h+(\delta r{)}^2+M\tau ⟨A_h{\mu }^{n+1},{\mu }^{n+1}{⟩}_h={\tilde{E}}_h(u^n,r^n),$ where $\delta u=u^{n+1}-u^n$, $\delta r=r^{n+1}-r^n$, and ${\tilde{E}}_h(u,r)=\frac{\kappa }{2}⟨A_{h}u,u{⟩}_h+r^2\mathrm{.}$ Thus ${\tilde{E}}_h$ is nonincreasing. This conclusion is about the modified energy; it does not assert monotonicity of the original physical energy ${\mathcal{E}}_h(u)=\frac{\kappa }{2}⟨A_{h}u,u{⟩}_h+E_{1,h}(u)$.

Use the finite-dimensional realization and the hypotheses of Backward-Euler convex splitting for Allen–Cahn and Cahn–Hilliard. Assume $P_{c,h}$ is convex and coercive on $X_h$, and $P_{e,h}$ is differentiable and convex. For every $u^n$, $M>0$, and $\tau >0$, both branches have a unique minimizer. Their first-order conditions are exactly $\frac{u^{n+1}-u^n}{\tau }=-M{\mu }^{n+1}\text{(Allen–Cahn)},$ $\frac{u^{n+1}-u^n}{\tau }=-M{A}_h{\mu }^{n+1}\text{(Cahn–Hilliard)}\mathrm{.}$ The Cahn–Hilliard branch conserves $⟨u^n,\mathbf{1}{⟩}_h$. With ${\mathcal{E}}_h$ as in the algorithm, ${\mathcal{E}}_h(u^{n+1})+M\tau \mathrm{\parallel }{\mu }^{n+1}{\mathrm{\parallel }}_h^2\le {\mathcal{E}}_h(u^n)\text{(Allen–Cahn)},$ ${\mathcal{E}}_h(u^{n+1})+M\tau ⟨A_h{\mu }^{n+1},{\mu }^{n+1}{⟩}_h\le {\mathcal{E}}_h(u^n)\text{(Cahn–Hilliard)}\mathrm{.}$ For the quartic split in Allen–Cahn and Cahn–Hilliard double-well gradient flows, the required convexity and coercivity hold for both realizations.