Abstract
We propose a robust and efficient augmented Lagrangian-type preconditioner for solving linearizations of the Oseen–Frank model arising in nematic and cholesteric liquid crystals. By applying the augmented Lagrangian method, the Schur complement of the director block can be better approximated by the weighted mass matrix of the Lagrange multiplier, at the cost of making the augmented director block harder to solve. In order to solve the augmented director block, we develop a robust multigrid algorithm which includes an additive Schwarz relaxation that captures a pointwise version of the kernel of the semi-definite term. Furthermore, we prove that the augmented Lagrangian term improves the discrete enforcement of the unit-length constraint. Numerical experiments verify the efficiency of the algorithm and its robustness with respect to problem-related parameters (Frank constants and cholesteric pitch) and the mesh size.
Similar content being viewed by others
1 Introduction
Liquid crystals (LC), first discovered by Reinitzer in 1888 [46], are materials that can exist in an intermediate mesophase between isotropic liquids and solid crystals: they can flow like liquids while also possessing long-range orientational order. Based on different ordering symmetries, Friedel [25] proposed to classify them into three broad categories: nematic, smectic and cholesteric. The nematic phase is the simplest and most extensively studied form of LC, where the molecules locally tend to align in one preferred direction, described in this work by a director field \(\mathbf {n}: \varOmega \rightarrow \mathbb {R}^3\). In the smectic phase, the molecules exhibit orientational order but also organize themselves into well-defined layers that can slide over each other. In the cholesteric phase, also referred as the chiral nematic phase, the molecules are arranged in layers, each of which is rotated with a fixed angle relative to the previous one. The distance over which the layers rotate by \(2\pi \) is referred to as the cholesteric pitch \(q_0\). A nonzero parameter \(q_0\) indicates chirality, while a zero value of \(q_0\) represents a nematic phase. Since the orientational properties of LC can be manipulated by imposing electric fields, they are often used to control light and have formed the basis of several important technologies in the area of display devices. Several thorough overviews on LC modeling and its history can be found in [5, 12, 52].
There are several models describing LC, e.g., the Oseen–Frank, Ericksen and Landau–de Gennes theories. The Oseen–Frank model [24, 42] is commonly used for the equilibrium orientation of liquid crystals. It employs a director \(\mathbf {n}: \varOmega \rightarrow \mathbb {R}^3\) as the state variable and minimizes a free energy functional. By definition, the director is a unit vector denoting the average orientation of the molecules in a fluid element at a point and headless in the sense that \(\mathbf {n}\) and \(-\mathbf {n}\) are indistinguishable. The free energy functional depends on Frank constants that describe the relative energetic costs of various kinds of distortions. We refer to [15, 19] for other continuum models such as the Ericksen and the Landau–de Gennes models. In this work, we will focus on the continuum Oseen–Frank theory. The key difficulty is that enforcing the unit-length constraint \(\mathbf {n}\cdot \mathbf {n} = 1\) with a Lagrange multiplier leads to a saddle-point system, which poses challenges because of its poor spectral properties. Several classical techniques regarding the solution of saddle-point problems are reviewed and illustrated in [10].
There are several existing works concerning preconditioners for Oseen–Frank models of nematic LC. For the saddle-point structure of harmonic maps (arising when all Frank constants are equal), Hu et al. [32] propose to use a block-diagonal preconditioner, consisting of a symmetric and spectrally equivalent multigrid operator and a discrete Laplacian operator. Ramage and Gartland [44] consider the case of an electrically coupled equal-constant nematic LC and combine a discretize-then-optimize approach with projection onto the nullspace of the discrete constraint to reduce the size of the linear system. The projected problem is then preconditioned with a block-diagonal preconditioner. Furthermore, a number of other preconditioners are discussed and analyzed in [7, 8] for double saddle-point systems arising in both potential fluid flows and electric-field coupled nematic LC. Concerning the double saddle-point structure, a class of Uzawa-type methods, which can be interpreted as generalized Gauss–Seidel methods, and an augmented Lagrangian technique are studied in [9]. It is shown that the applied augmented Lagrangian form is mesh-independent and the performance of the iteration can be improved by increasing the value of \(\gamma \). These references also apply the discretize-then-optimize approach to tackle the pointwise unit-length vector constraint. In this paper, we will employ the optimize-then-discretize strategy and enforce the unit-length constraint on the continuous level. As an alternative to block preconditioning strategies, monolithic multigrid methods for the nematic problem have been proposed using Vanka [1] and Braess–Sarazin [3] relaxation.
There is less work on preconditioning for cholesteric LC. A damped Newton method with LU decomposition was applied to the bifurcation analysis of cholesteric problem in [18] with good results, but no discussion of preconditioners is presented.
In this paper, we propose to enforce the unit-length constraint with an augmented Lagrangian approach to help control the Schur complement arising in the saddle-point system. When combined with specialized multigrid schemes, augmented Lagrangian strategies can yield scalable, mesh-independent, and parameter-robust preconditioners. A notable success is the development of Reynolds-robust solvers for the two- [11, 41] and three-dimensional [21] stationary Navier–Stokes equations.
This success motivates the investigation of whether similar ideas can underpin robust solvers in the LC case.
The main contribution of this work is the development of a robust multigrid solver for the augmented director block and an effective Schur complement approximation for the linearization of the cholesteric Oseen–Frank equations. The robust multigrid strategy is motivated by the general theory of Schöberl and Lee et al. [35, 49, 50]. We develop a multigrid relaxation scheme that captures an approximation to the kernel of the semi-definite augmentation term and account for this approximation in the spectral analysis. Furthermore, a proof of the improvement of the discrete constraint is given and verified numerically. A key difference to previous applications of these ideas in linear elasticity and the Navier–Stokes equations is that the constraint to be imposed on the director is nonlinear.
This paper is organized as follows. The Oseen–Frank model is reviewed in Sect. 2 and the solvability of the discretized Newton linearizations is briefly analyzed. The augmented Lagrangian strategy for enforcing the unit-length constraint is discussed. A Picard iteration is proposed for solving the augmented nonlinear equations. We then give a theoretical justification of the continuous and discrete augmented Lagrangian stabilizations in Sect. 3. This further leads to our choice of the approximation to the Schur complement matrix arising from the Picard iteration. In Sect. 4, we prove that the augmented Lagrangian strategy improves the discrete enforcement of the constraint. A robust multigrid algorithm for the augmented top-left block is discussed in Sect. 5 which also includes a formal spectral analysis of our preconditioner with the property of the approximate kernel. Numerical experiments in two-dimensional domains are reported in Sect. 6 to verify the effectiveness and robustness of our proposed augmented Lagrangian preconditioner. Finally, some conclusions are presented in Sect. 7.
2 Oseen–Frank model
Let \(\varOmega \subset \mathbb {R}^d, d=\{2,3\}\) be an open, bounded domain with Lipschitz boundary \(\partial \varOmega \) and denote \(\mathbf {H}^1_g(\varOmega ) = \{\mathbf {v} \in H^1(\varOmega ; \mathbb {R}^3): \mathbf {v}|_{\partial \varOmega } = \mathbf {g}\}\) with a vector field \(\mathbf {g}\in H^{1/2}(\partial \varOmega ;\mathbb {S}^2)\). Here, \(\mathbb {S}^2\) represents the surface of the unit ball centered at the origin. Assume that the cholesteric LC occupying the domain \(\varOmega \) is equipped with a rigid anchoring (Dirichlet) boundary condition \(\mathbf {n}|_{\partial \varOmega }=\mathbf {g}\)Footnote 1. The Oseen–Frank model (cf. [24]) considers the following minimization problem:
where the Frank energy density \(W(\mathbf {n})\) is of the form
where \(\text {tr}(\cdot )\) denotes the trace of a matrix, \(K_i \in \mathbb {R}\) \((i=1,2,3,4)\) are elastic constants (called Frank constants) and \(q_0 \ge 0\) is the preferred pitch for the cholesteric. \(K_1\), \(K_2\), \(K_3\), and \(K_4\) are referred to as the splay, twist, bend, and saddle-splay constants, respectively. Note here \(\nabla \mathbf {n}\) is matrix-valued and \((\nabla \mathbf {n})^2\) denotes the matrix multiplication of the matrix \(\nabla \mathbf {n}\) and itself.
If \(K_1=K_2=K_3=K>0\) and \(K_4=0\), the energy density (2.2) reduces to the so-called equal-constant approximation, with energy density
which is a useful simplification to help us gain qualitative insight into more complex situations.
Remark 2.1
When \(q_0=0\), the energy density (2.2) corresponds to the nematic case. Furthermore, when combined with the equal-constant approximation, (2.2) reduces to
With this free energy density, the solution to the minimization problem (2.1) is unique and is known as the harmonic map from a two- or three-dimensional compact manifold to \(\mathbb {S}^2\) [36]. Some fast numerical algorithms for (2.3) have been proposed and tested in [32].
The last term (the saddle-splay term or the null Lagrangian) in (2.2) can be dropped as its integral reduces to a surface integral, which is essentially a constant if applying Dirichlet boundary conditions to the model, via the divergence theorem. For mixed periodic and Dirichlet boundary conditions considered in Sect. 6.2.1, we can verify directly that this saddle-splay energy vanishes. Hence, for simplicity, it suffices to consider the following Frank energy density
In this paper, we use a more compact form of the free energy (2.1) as in [2, 3] by introducing a symmetric dimensionless tensor
where \(\kappa = K_2/ K_3\) and \(\mathbf {I}\) is the second-order identity tensor. By the classical equality
the original energy functional \(J(\mathbf {n})\) can be written as
Here and throughout this work, \(\langle \cdot , \cdot \rangle _0\) denotes the inner product in \(L^2(\varOmega )\) with its induced norm \(\Vert \cdot \Vert _0\). It can be observed that the auxiliary tensor \(\mathbf {Z}\) contributes to the nonlinearity of \(J(\mathbf {n})\) in (2.5).
Remark 2.2
There is another widely used simplification of the energy density (2.2), where \(q_0=0\) and \(K_2=K_3=K_1+K\), \(K_4=-K\) [29, 37]. In this case, (2.2) becomes
and it is expected that as \(K\rightarrow \infty \), the asymptotic behavior of minimizers provides a description of the phase transition process of LC from the nematic to the smectic-A phases [29, 37, 38].
Furthermore, it is proven in [2, Section 2.3] that \(\mathbf {Z}\) is uniformly (with respect to \(\mathbf {x}\in \varOmega \)) symmetric positive definite (USPD) as long as sufficient control is maintained on \(\mathbf {n}\cdot \mathbf {n}-1\). This property of \(\mathbf {Z}\) plays an essential role in proving the well-posedness of the saddle-point problem in the nematic case. We restate the result of \(\mathbf {Z}\) being USPD in the following, as it is important later:
Lemma 2.1
[2, Section 2.3] Assume \(\alpha \le |\mathbf {n}|^2 \le \beta \) \(\forall \mathbf {x}\in \varOmega \) with \(0<\alpha \le 1\le \beta \). If \(\kappa >1\), then \(\mathbf {Z}\) is USPD on \(\varOmega \); for \(0<\kappa <1\), then \(\mathbf {Z}\) is USPD on \(\varOmega \) if \(\beta <\frac{1}{1-\kappa }\).
Remark 2.3
Notice that the regularity of \(\mathbf {n}\in \mathbf {H}^1(\varOmega )\) is enough for the functional \(J(\mathbf {n})\) of (2.5) to be well defined. In fact, \(\mathbf {n}\in \mathbf {H}^1(\varOmega )\) implies \(\nabla \cdot \mathbf {n}\) and \(\nabla \times \mathbf {n}\) in \(\mathbf {L}^2(\varOmega )\). By (2.4), \(\mathbf {n}\cdot (\nabla \times \mathbf {n})\in L^2(\varOmega )\). This ensures that the term \(\langle q_0, \mathbf {n}\cdot (\nabla \times \mathbf {n})\rangle _0\) in (2.5) is defined. Furthermore, Lemma 2.1 gives the boundedness of \(\mathbf {Z}\), which guarantees the \(L^2\)-regularity of the term \(\mathbf {Z}\nabla \times \mathbf {n}\) in (2.5).
Naturally, the values of elastic constants and the cholesteric pitch will be an important factor in determining the minimizers. In order to satisfy non-negativity of the energy density, i.e.,
we need additional assumptions on those constants. This gives rise to Ericksen’s inequalities (see [5, 6] and references therein):
Remark 2.4
We have included the inequalities with regard to constant \(K_4\) here for generality, though they are not necessary in our work as we have eliminated the \(K_4\)-related term in the free energy. In this paper, we will simply consider \(K_i>0\) (\(i=1,2,3\)) to avoid any technical issues.
For the minimization problem (2.1) arising in (nematic or cholesteric) liquid crystals, it has been proven in [36, Theorem 2.1] that there exists a solution.
Theorem 2.1
[36, Theorem 2.1] Let \(\varOmega \) be a bounded Lipschitz domain and assume the Dirichlet boundary data \(\mathbf {g}\in {H}^{1/2}(\partial \varOmega ;\mathbb {S}^2)\). If \(K_1,K_2,K_3>0\), then there exists an \(\mathbf {n}\in {H}_g^1(\varOmega ;\mathbb {S}^2):=\{\mathbf {n}\in H^1(\varOmega ;\mathbb {S}^2):\mathbf {n}=\mathbf {g}\ \text {on }\partial \varOmega \}\) such that
The main difficulty in solving the Oseen–Frank model (2.1) is the enforcement of the unit-length constraint. There are several existing approaches to handling constraints, e.g., projection [38], Lagrange multipliers, and penalty methods [39, Section 12.3 & 17].
The projection method is numerically simple but the value of the energy functional may go up and down dramatically after each projection, making it difficult to control in the optimization procedure [38]. A Lagrange multiplier is often used to replace constrained optimization problems with unconstrained ones, but an important disadvantage of this approach is that it introduces another unknown (i.e., the Lagrange multiplier) and leads to a saddle-point structure which can be difficult to solve [10]. On the other hand, the penalty method has the favorable property that the resulting system has an energy decay property [37] which may result in an easier theoretical and numerical study of the solution. However, the penalty parameter has to be very large for the accuracy of approximating the constraints, leading to an ill-conditioned system. Some works based on either projection or pure penalty methods for nematic phases can be found in [28, 29, 37] and the references therein.
Fortunately, it is possible to amend the ill-conditioning effects with large penalty parameters that are inherent in the pure penalty method by combining it with a Lagrange multiplier. This is the augmented Lagrangian (AL) algorithm [23]. This strategy combines the advantages of both schemes: the penalty parameter can be relatively small due to the presence of the Lagrange multiplier, and the Schur complement of the saddle-point system is easier to solve due to the presence of the penalty term [11, 21, 28, 29, 40].
We first consider the method of Lagrange multipliers. We then add the augmented Lagrangian term to control the Schur complement of the system.
2.1 Lagrange multiplier and Newton linearization
By introducing the Lagrange multiplier \(\lambda \in L^2(\varOmega )\), the associated Lagrangian of the minimization problem (2.1) is then defined as
and its first-order optimality conditions are: find \((\mathbf {n}, \lambda ) \in \mathbf {H}^1_g(\varOmega ) \times L^2(\varOmega )\) such that
As (2.7) is nonlinear, Newton linearization is employed. Let \(\mathbf {n}_k\) and \(\lambda _k\) be the current approximations for \(\mathbf {n}\) and \(\lambda \), respectively, and denote the corresponding updates to these approximations as \(\delta \mathbf {n} = \mathbf {n}_{k+1} - \mathbf {n}_k\) and \(\delta \lambda = \lambda _{k+1} - \lambda _k\). Then the Newton iteration at \((\mathbf {n}_k,\lambda _k)\) in block form is given by: find \((\delta \mathbf {n},\delta \lambda )\in \mathbf {H}^1_0(\varOmega )\times L^2(\varOmega )\) such that
where
and
Since \(\mathscr {L}(\mathbf {n}, \lambda )\) is linear in \(\lambda \), \(\mathscr {L}_{\lambda \lambda } = 0\). This results in (2.8) being a saddle-point problem.
With a suitable spatial discretization (we only consider conforming finite elements in this work, i.e., \(V_h\subset \mathbf {H}^1_0(\varOmega )\), \(Q_h\subset L^2(\varOmega )\)), a symmetric saddle-point system must be solved at each Newton iteration:
where U and P represent the coefficient vectors of \(\delta \mathbf {n}\) and \(\delta \lambda \) in terms of the basis functions of \(V_h\) and \(Q_h\), respectively.
We can accordingly write the discrete variational problem as: find \(\delta \mathbf {n}_h\in V_h\) and \(\delta \lambda _h\in Q_h\) such that
where \(a(\cdot , \cdot )\) and \(b(\cdot ,\cdot )\) are bilinear forms given by
and
and F and G are linear functionals in the forms of
and
Remark 2.5
The well-posedness of the continuous and discretized Newton system (with the \(([\mathbb {Q}_m]^d\oplus B_F)\)-\(\mathbb {Q}_0\) finite element pair, \(m\ge 1\)) for a generalized nematic LC problem is discussed in [2], where \(B_F\) denotes the space of quadratic bubbles and \(\mathbb {Q}_k\) represents tensor product piecewise \(C^0\) polynomials of degree \(k \ge 0\) on a quadrilateral mesh. Moreover, the authors of [3] considered the pure penalty approach for nematic LC and obtained a well-posedness result of the penalized Newton iteration through similar techniques. We will follow these analysis strategies in this section.
In our work, we will denote by \(\mathbb {P}_k\) the set of piecewise \(C^0\) polynomials of degree \(k\ge 0\) on a mesh of triangles or tetrahedra.
It is straightforward to deduce the well-posedness of the discrete Newton iteration (2.11) for cholesteric problems under some proper assumptions on the problem-dependent constants. In fact, two additional \(q_0\)-related terms in \(\mathscr {L}_{\mathbf {n}\mathbf {n}}\) from (2.9) compared to the nematic energy density from [2] are simply \(L^2\) inner products, which can be easily bounded using the Cauchy–Schwarz and triangle inequalities. We state the results without proof in the following and start with some assumptions.
Assumption 2.1
Assume that there exist constants \(0<\alpha \le 1\le \beta \) such that \(\alpha \le |\mathbf {n}_k|^2 \le \beta \). For \(0<\kappa <1\), assume further that \(\beta < \frac{1}{1-\kappa }\). By Lemma 2.1, \(\mathbf {Z}(\mathbf {n}_k)\) remains USPD with lower bound \(\eta \) and upper bound \(\Lambda \), i.e.,
Note that here and hereafter, \(\Vert \cdot \Vert _1\) denotes the \(H^1\) norm: \(\Vert w\Vert _1^2=\Vert w\Vert _0^2+\Vert \nabla w\Vert _0^2\).
Lemma 2.2
(Continuous coercivity) With Assumption 2.1, we assume further that the current Lagrange multiplier approximation \(\lambda _k\) is pointwise non-negative almost everywhere. Let \(K_1>K_2q_0C_4\) and \(K_3\eta > K_2q_0(C_4 + 1)\) with \(C_4\) to be defined. Then there exists an \(\alpha _0>0\) such that
Moreover, when \(\kappa =1\), i.e., \(K_2=K_3\), if \(K_1>K_2q_0C_4\) and \(1>q_0(C_4 + 1)\), then the coercivity result (2.12) also holds.
Proof
With the lower bound \(\eta \) of \(\mathbf {Z}\), we compute the bilinear form:
where the first inequality comes from the assumption that \(\lambda _k\) is non-negative pointwise and the last two inequalities are derived by Cauchy–Schwarz and Hölder inequalities, respectively.
By Remark 2.7 of [27], for a bounded Lipschitz domain, there exists \(C_1>0\) such that
for all \(\mathbf {v}\in \mathbf {H}_0(\mathrm {div},\varOmega )\cap \mathbf {H}_0(\mathrm {curl},\varOmega )\).Footnote 2 Here, we denote
where \(\varvec{\nu }\) is the outward unit normal on the boundary \(\partial \varOmega \). Then using the classical Poincaré inequality, \(\Vert \mathbf {v}\Vert _0^2\le C_3\Vert \nabla \mathbf {v}\Vert _0^2\) for all \(\mathbf {v}\in \mathbf {H}^1_0(\varOmega )\), and defining \(C_4=C_1C_3>0\), we have
Furthermore, there exists \(C_2=C_4+C_1>0\) such that
It follows that
Choosing \(c=\min \{K_1-K_2q_0C_4, K_3\eta -K_2q_0C_4-K_2q_0 \}>0\) (the positivity follows from the assumptions) and \(\alpha _0=c/C_2\), we find that the coercivity (2.12) holds.
In particular, when \(\kappa =1\) (i.e., \(K_2=K_3\)), we have \(\mathbf {Z}=\mathbf {I}\) and thus \(\eta =1\). Then, the bilinear form becomes
By choosing \(C = \min \{K_1-K_2q_0 C_4, K_2(1 - q_0 C_4 - q_0)\}>0\) (the positivity comes from the assumptions) and \(\alpha _0=C/C_2\), we obtain the desired coercivity
as stated in (2.12). \(\square \)
So far, the coercivity of the bilinear form \(a(\cdot ,\cdot )\) has been shown for all functions in \(\mathbf {H}^1_0(\varOmega )\). The discrete coercivity follows if a conforming finite element for the director space is chosen.
The boundedness of the bilinear form \(a(\cdot ,\cdot )\) and the right hand side functionals \(F(\cdot )\) and \(G(\cdot )\) can be obtained directly by following the proofs in [2]. Hence, we omit the details here.
It remains to consider the discrete inf-sup condition of the bilinear form \(b(\cdot ,\cdot )\) for a finite element pair \(V_h\)-\(Q_h\), i.e. whether there exists a constant c such that
The continuous inf-sup condition was shown in [17, Appendix B] and [32, Theorem 3.1]. However, the discrete inf-sup condition is not inherited from the continuous problem. Some previous works have succeeded in obtaining a discrete inf-sup condition for some specific discretizations. A discrete inf-sup condition was proven for the \(([\mathbb {Q}_m]^d\oplus B_F)\)-\(\mathbb {Q}_0\) element on quadrilaterals in [17, Lemma 2.5.14] and [2, Lemma 3.12]. The discrete inf-sup condition for the \([\mathbb {P}_1]^2\)-\(\mathbb {P}_1\) discretization is shown in [32, Theorem 4.5], where the analysis is only valid for the two-dimensional case due to the use of some special inverse inequalities. It is straightforward to deduce that an enrichment of \(V_h\) still guarantees the stability of the discretization, and thus \([\mathbb {P}_2]^2\)-\(\mathbb {P}_1\) is inf-sup stable under the same conditions.
We now consider the matrix form of the saddle-point system (2.10). The coercivity of the bilinear form \(a(\cdot ,\cdot )\) implies the invertibility of the coefficient matrix A and the discrete inf-sup condition indicates that B has full row rank. We use the full block factorization preconditioner
with approximate inner solves \(\tilde{A}^{-1}\) and \(\tilde{S}^{-1}\) for the director block and the Schur complement \(S = - BA^{-1}B^\top \), respectively, for solving the saddle-point problem (2.10). With exact inner solves, this is an exact inverse. With this strategy, solving the original saddle-point problem (2.10) reduces to solving two smaller linear systems involving A and S. Even though A is sparse, its inverse is generally dense, making it impractical to store S explicitly. In this situation, developing a fast solver for A is tractable while approximating S becomes difficult. We will return to this issue in Sects. 3 and 5.
2.2 Augmented Lagrangian form
Now, we employ the AL stabilization strategy and modify the linearized saddle point system to control its Schur complement S.
2.2.1 Penalizing the constraint
We penalize the continuous form of the nonlinear constraint \(\mathbf {n}\cdot \mathbf {n}=1\) in the AL algorithm and obtain the Lagrangian
for \(\gamma \ge 0\). The weak form of the associated first-order optimality conditions is to find \((\mathbf {n}, \lambda )\in \mathbf {H}^1_g(\varOmega )\times L^2(\varOmega )\) such that
The Newton linearization at a given approximation \((\mathbf {n}_k,\lambda _k)\) yields a system of the form:
Thus, we have to solve the augmented discrete variational problem:
where
and
Comparing (2.14) to the original system (2.11), only the bilinear form \(a(\cdot ,\cdot )\) and the right hand side functional \(F(\cdot )\) have changed. The boundedness of \(F^c(\cdot )\) follows straightforwardly via the Cauchy–Schwarz inequality. As for the coercivity of \(a^c(\cdot ,\cdot )\), an additional assumption on the penalty parameter \(\gamma \) is needed.
Lemma 2.3
(Continuous coercivity) Let \(\alpha _0 > 0\) be the coercivity constant of \(a(\cdot , \cdot )\). If \(\alpha _0>2\gamma |\alpha -1|\) with \(0<\alpha \le 1\le \beta \) satisfying \(\alpha \le |\mathbf {n}_k|^2\le \beta \), there exists a \(\beta _0>0\) such that
Proof
Note that
By the assumption that \(a(\mathbf {v},\mathbf {v})\ge \alpha _0\Vert \mathbf {v}\Vert _1^2\) for some \(\alpha _0>0\), we have
Moreover, since \(\mathbf {n}_k\cdot \mathbf {n}_k\ge \alpha \) and \(\alpha - 1\le 0\), we get
Thus, by taking \(\beta _0=\alpha _0-2\gamma |\alpha -1|>0\), we obtain the desired coercivity property. \(\square \)
The condition \(\alpha _0 > 2\gamma |\alpha -1|\) in Lemma 2.3 indicates a limit on the value of \(\gamma \) to ensure the solvability of the augmented system (2.14). However, it is desirable to use large values of \(\gamma \) to achieve better control of the Schur complement. We therefore choose to employ a Picard iteration to solve the nonlinear problem, omitting the term \(2\gamma \langle \mathbf {n}_k\cdot \mathbf {n}_k-1, \mathbf {v}\cdot \mathbf {v}\rangle _0\) from the linearized equations. This yields the linearized problem: find \((\delta \mathbf {n}_h, \delta \lambda _h) \in V_h \times Q_h\) such that
with the modified bilinear form
to be solved at each nonlinear iteration. This ensures that the (1, 1)-block is coercive with a coercivity constant independent of \(\gamma \). Moreover, in contrast to the situation with the Navier–Stokes equations, numerical experiments indicate that the use of this Picard requires fewer nonlinear iterations to converge to a given tolerance than using the full Newton linearization (see Sect. 6.2.1).
The corresponding matrix form of the variational problem (2.14) becomes
where \(A_*\) is the assembly of \(4\langle \mathbf {n}_k\cdot \mathbf {u}, \mathbf {n}_k\cdot \mathbf {v}\rangle _0\) and l denotes the assembly of \(- 2\langle \mathbf {n}_k\cdot \mathbf {n}_k-1, \mathbf {n}_k\cdot \mathbf {v}\rangle _0\). Note that compared to the non-augmented version (2.10), the (1, 1) block in (2.17) has an additional semi-definite term \(\gamma A_*\) with a large coefficient \(\gamma \). Its sparsity pattern remains unchanged. We will construct a robust multigrid method to solve this top-left block in Sect. 5.
Since the unit-length constraint is enforced exactly in (2.13), the continuous solutions to minimizing both (2.13) and (2.6) are the same. However, the unit-length constraint is not enforced exactly in our finite element discretization, and hence this stabilization does change the computed discrete solution.
Remark 2.6
When applying the augmented Lagrangian strategy, one can apply it before discretization or afterwards. In this work we apply the continuous penalization, as it improves the enforcement of the nonlinear constraint, as shown later in Sect. 4. This is different to the approach considered in [11, 21] for the stationary Navier–Stokes equations, where the discrete AL stabilization was used to yield a system that has the same solution but a better Schur complement.
3 Approximation to the Schur complement
The Schur complement of the augmented director block in (2.17) is given by
We now proceed to analyze this Schur complement by following similar techniques to those of [31, §4]. We will show that \(A_*\) is equal to the matrix arising from the discrete AL stabilization (which controls the Schur complement) plus a perturbation that vanishes as the mesh is refined.
Let \(\varPi _{Q_h}: Q\rightarrow Q_h\) be the orthogonal \(L^2\) projection operator, i.e.,
We define the fluctuation operator \(\kappa :={I}-\varPi _{Q_h}\) where \({I}:Q\rightarrow Q\) is the identity mapping. Therefore, one has
For \(\mathbf {u}_h,\mathbf {v}_h\in V_h\), one can split the term \(4\langle \mathbf {n}_k \cdot \mathbf {u}_h, \mathbf {n}_k \cdot \mathbf {v} \rangle _0\) into the following terms using the properties of \(\kappa \) and \(\varPi _{Q_h}\):
Note here that the assembly of the first term is \(B^\top M_\lambda ^{-1}B\), where \(M_\lambda \) is the mass matrix associated with the finite element space for the multiplier \(Q_h\). This can then be readily used with the Sherman–Morrison–Woodbury formula to derive an approximation of the Schur complement. The second term \(\langle \kappa (2\mathbf {n}_k\cdot \mathbf {u}),\kappa (2\mathbf {n}_k\cdot \mathbf {v})\rangle _0\) characterizes the difference between \({A}_*\) and \(B^\top M_\lambda ^{-1}B\). The next result shows that it vanishes as the mesh size \(h\rightarrow 0\) (see Theorem 3.1) and thus, in this limit, the tractable term \(B^\top M_\lambda ^{-1}B\) dominates \({A}_*\).
Theorem 3.1
Let \((\delta \mathbf {n}_h, \delta \lambda _h)\in V_h\times Q_h\) be the solution of the augmented discrete system (2.15) with corresponding degrees of freedom \((U, P)\in \mathbb {R}^n\times \mathbb {R}^m\). Then, for the Newton linearization at a given approximation \((\mathbf {n}_k, \lambda _k)\) satisfying \(\alpha \le |\mathbf {n}_k|^2\le \beta \) with \(0< \alpha \le 1\le \beta \) and \(|\nabla \mathbf {n}_k|\) bounded pointwise a.e., we have
where \(\Vert \cdot \Vert _{\mathbb {R}^n}\) denotes the Euclidean norm.
Proof
Assuming \(\mathbf {v}_h\in V_h\) and using the basis representations in \(V_h=\mathrm {span}\{\phi _i\}\) for \(\delta \mathbf {n}_h\) and \(\mathbf {v}_h\):
we obtain
by applying the Cauchy–Schwarz inequality.
One readily sees that \(G_1\le C_1\) for a certain constant \(C_1\) from the continuity of \(\kappa \). Furthermore, we write
Note that [34, Theorem 3.43] as used in [31] gives the relation between the discrete vector Y and its associated continuous function \(\mathbf {v}_h\):
for some \(C_r>0\). Then with the fact that \(\mathbf {n}_k\) is bounded we have
Moreover, [14, Theorem 1] implies
we deduce the following \(L^2\)-projection error estimate
Note here we have used the pointwise boundedness of \(\mathbf {n}_k, \nabla \mathbf {n}_k\) a.e. and the fact that \(\delta \mathbf {n}_h\in V_h\subset H^1(\varOmega )\).
Combining these estimates regarding \(G_1,G_2, G_3\), we find
\(\square \)
This result suggests the use of the algebraic approximation
The reason for doing so is that we can straightforwardly calculate the inverse of this approximation (3.1) by the Sherman–Morrison–Woodbury formula as follows:
The solver requires the action of \(S_\gamma ^{-1}\), i.e., solving linear systems involving \(S_\gamma \). For large \(\gamma \), a simple and effective approach is to employ the approximation
On the infinite-dimensional level, the effect of the augmented Lagrangian term is to make \(-\gamma ^{-1} I\) (I the identity operator on the multiplier space) an effective approximation for the Schur complement [43, Lemma 3]. When discretized, this indicates that the weighted multiplier mass matrix \(-\gamma ^{-1} M_\lambda \) will be an effective approximation for \(S_\gamma \), with the approximation improving as \(\gamma \rightarrow \infty \).
In fact, the approximation of the inverse of the discretely augmented Schur complement (3.2) can be improved further by combining \(-\gamma M_{\lambda }^{-1}\) with a good approximation of the unaugmented Schur complement S [30]. Given an approximation \(\tilde{S}\) of S, we employ
It is therefore of interest to consider the Schur complement of the unaugmented problem. In the context of the Stokes equations, the Schur complement is spectrally equivalent to the viscosity-weighted pressure mass matrix [16, 51, 54]. Following similar techniques, an approximation can be obtained by proving that \(BA^{-1}B^\top \) is spectrally equivalent to \(M_{\lambda }\) for the equal-constant nematic case. This gives us good insight into the choice of \(\tilde{S}^{-1}\).
Theorem 3.2
For equal-constant nematic LC problems without augmented Lagrangian stabilization, the matrix \(BA^{-1}B^\top \) arising from the Newton-linearized system is spectrally equivalent to the multiplier mass matrix \(M_\lambda \), under the same assumptions as in Lemma 2.2.
Proof
For the equal-constant model with Dirichlet boundary conditions \(\mathbf {n}=\mathbf {g}\in H^{1/2}(\partial \varOmega ;\mathbb {S}^2)\), its corresponding Lagrangian is
After Newton linearization and introducing conforming finite dimensional spaces \(V_h\subset \mathbf {H}^1_0(\varOmega )\) and \(Q_h \subset L^2(\varOmega )\), the discrete variational problem is to find \(\delta \mathbf {n}_h\in V_h\), \(\delta \lambda _h\in Q_h\) satisfying
where \(\mathbf {n}_k\) and \(\lambda _k\) represent the current approximations to \(\mathbf {n}\) and \(\lambda \), respectively. This can be rewritten in block matrix form as
where as before \(U\in \mathbb {R}^n\) and \(P\in \mathbb {R}^m\) are the unknown coefficients of the discrete director update and the discrete Lagrange multiplier update with respect to the basis functions in \(V_h\) and \(Q_h\), and A denotes the symmetric form \(K\langle \nabla \delta \mathbf {n}_h,\nabla \mathbf {v}_h\rangle _0 + 2\langle \lambda _k,\delta \mathbf {n}_h\cdot \mathbf {v}_h\rangle _0\). The coercivity property of the bilinear form from Lemma 2.2 ensures that A is positive definite.
The coefficient matrix \(\mathscr {A}\) is symmetric and indefinite (resulting in \(\mathscr {A}\) possessing both positive and negative eigenvalues). Moreover, \(\mathscr {A}\) is non-singular if and only if B has full row rank, which can be deduced from the discrete inf-sup condition.
Denote
Notice that the validity of the first norm follows from the assumed pointwise non-negativity of \(\lambda _k\).
For a stable mixed finite element, from the inf-sup condition, there exists a positive constant c independent of the mesh size h such that
leading to its matrix form
Thus, we have
where the maximum is attained at \(z=(P^\top BA^{-1/2})^\top \). It yields
Regardless of the stability of the finite element pair, we can deduce from the boundedness of B that there exists a positive constant \(c_1\) such that
Hence,
where again the maximum is attained at \(z=(P^\top BA^{-1/2})^\top \). This gives rise to
Therefore for inf-sup stable finite element pairs, we have by (3.4) and (3.5)
This indicates that \(BA^{-1}B^\top \) is spectrally equivalent to \(M_\lambda \). \(\square \)
Remark 3.1
It follows from Theorem 3.2 that \(\gamma =0\) should show mesh-independence (i.e., the average number of FGMRES iterations per Newton iteration does not deteriorate as one refines the mesh) in the case of equal-constant nematic LC. This can be observed in subsequent numerical experiments reported in Table 6 (see the column where \(\gamma =0\)). One should also notice that such mesh-independence for \(\gamma =0\) is also shown in Table 2 for the non-equal-constant case, suggesting it has use outside the context of augmented Lagrangian methods also.
Combining Theorem 3.2 with (3.3), our final approximation for \(S_\gamma ^{-1}\) is given by
4 Improvement of the constraint
We have now observed that the continuous AL form introduced in Sect. 2.2.1 can help control the Schur complement. Another contribution of this AL stabilization is that it improves the discrete constraint as we increase the value of the penalty parameter \(\gamma \). An example of improving the linear divergence-free constraint in the Stokes system can be found in [33, Section 5.1]. In this section, we will use a similar strategy to show the improvement of the discrete constraint as \(\gamma \) increases.
We restrict ourselves to the equal-constant case with constant Dirichlet boundary conditions. That is to say, we consider the Oseen–Frank model with Dirichlet boundary condition \(\mathbf {n}|_{\partial \varOmega }=\mathbf {g}\), where \(\mathbf {g}\) is a nonzero constant vector satisfying \(|\mathbf {g}|=1\). We use the \([\mathbb {P}_1]^d\)-\(\mathbb {P}_1\) finite element pair in this section, so both the director \(\mathbf {n}\) and the Lagrange multiplier \(\lambda \) are approximated by continuous piecewise-linear polynomials. For this section, we denote finite element spaces for the director and the Lagrange multiplier by \(V_{h,g}:=V_h \cap \mathbf {H}^1_g(\varOmega )\) and \(Q_h \subset L^2(\varOmega )\), respectively, and denote \(V_{h,0}= V_h\cap \mathbf {H}^1_0(\varOmega )\).
We restate the associated nonlinear discrete variational problem as follows: find \((\mathbf {n}_h, \lambda _h)\in V_{h,g}\times Q_h\) such that
Take the test function \(\mathbf {v}_h = \mathbf {n}_h-{\mathbf {g}}\in V_{h,0}\) in (4.1a) to obtain
Note that in this step we have used the fact that since \(\mathbf {g}\) is a constant vector, its derivative is zero.
As (4.1b) is valid for arbitrary \(\mu _h\in Q_h\) and one can easily verify that \(\mathbf {n}_h\cdot \mathbf {g}\in Q_h\), we have
Then taking \(\mu _h=1\) and \(\mu _h=\lambda _h\) leads to
respectively. Thus, (4.2) collapses to
By the Cauchy–Schwarz and Hölder inequalities, we observe an upper bound for the right hand side of (4.3):
Meanwhile, the left hand side of (4.3) can be bounded from below:
where \(|\varOmega |\) denotes the measure of the domain \(\varOmega \).
Hence, by combining (4.4) and (4.5), we have
Since the right hand side of (4.6) is a fixed constant independent of \(\gamma \), taking \(\gamma \) larger value forces the constraint approximation error \(\Vert \mathbf {n}_h\cdot \mathbf {n}_h-1\Vert _0\) to become smaller. In fact, (4.6) implies that \(\Vert \mathbf {n}_h\cdot \mathbf {n}_h-1\Vert _0\le \mathscr {O}(\gamma ^{-1/2})\).
Remark 4.1
The technique shown in this section can be extended in a similar way to the multi-constant case; we omit the details here for brevity.
5 A robust multigrid method for \(A_{\gamma }\)
As discussed in Sect. 3, the addition of the augmented Lagrangian term has the effect of controlling the Schur complement of the matrix in (2.17). However, the tradeoff is that it complicates the solution of the top-left block \(A_\gamma \), as it adds a semi-definite term with a large coefficient. For the augmented Lagrangian strategy to be successful, we require a \(\gamma \)-robust solver for the top-left block. Fortunately, a rich literature is available to guide the development of multigrid solvers for nearly singular systems [35, 49, 50]. In this section we develop a parameter-robust multigrid method for \(A_\gamma \).
Schöberl’s seminal paper on the construction of parameter-robust multigrid schemes [49] lists two requirements that must be satisfied for robustness. The first requirement is a parameter-robust relaxation method; this is achieved by developing a space decomposition that stably captures the kernel of the semi-definite terms. The second requirement is a parameter-robust prolongation operator, i.e. one whose continuity constant is independent of the parameters. This is achieved by (approximately) mapping kernel functions on coarse grids to kernel functions on fine grids. We discuss both of these requirements below.
For ease of notation, we consider the two-grid method applied to the equal-constant nematic case, and use subscripts h and H to distinguish fine and coarse levels respectively. That is to say, \(V_H\) represents the coarse-grid function space and \(A_{H,\gamma }:V_H\rightarrow V_H^*\) corresponds to the partial differential equations (PDEs) on \(V_H\).
For the domain \(\varOmega \), we consider a non-overlapping triangulation \(\mathscr {M}_H\), i.e.,
The fine grid \(\mathscr {M}_h\) with \(h=H/2\) is obtained by a regular refinement of the simplices in \(\mathscr {M}_H\). In what follows we consider both the \([\mathbb {P}_1]^d\)-\(\mathbb {P}_1\) and \([\mathbb {P}_2]^d\)-\(\mathbb {P}_1\) discretizations.
5.1 Relaxation
After applying the AL method introduced in Sect. 2.2.1, the discrete linear variational form corresponding to the top-left block \(A_{\gamma } = A+\gamma A_*\) is given by
with \(\mathbf {u}_h\in V_h\subset \mathbf {H}^1_0(\varOmega )\) being the trial function and \(\mathbf {v}_h\in V_h\) the test function. Note that \(\mathbf {n}_k\) and \(\lambda _k\) are the current approximations to the director \(\mathbf {n}\) and the Lagrange multiplier \(\lambda \), respectively, in the Newton iteration. The first two terms of \(a^m\) are symmetric and coercive because of the uniform non-negativity of \(\lambda _k\) in the assumption of our well-posedness result. The kernel of the semi-definite term involving \(\gamma \) is
In the case of \(\gamma \) being very large, the variational problem involving (5.1) is nearly singular and common relaxation methods like Jacobi and Gauss–Seidel will not yield effective multigrid cycles, as we explain below.
Relaxation schemes can be devised in a generic way by considering space decompositions
where the sum of vector spaces on the right is not necessarily a direct sum [56]. This space decomposition induces a relaxation method by (approximately) solving the Galerkin projection of the error equation onto each subspace \(V_i\), and combining the resulting estimates of the error. This can be done in an additive or multiplicative way. For example, if \(V_h = \mathrm {span}(\phi _1, \dots , \phi _N)\), Jacobi and Gauss–Seidel are induced by the space decomposition
where the updates are performed additively for Jacobi and multiplicatively for Gauss–Seidel. One of the key insights of [35, 49] was that the key requirement for parameter-robustness when applied to nearly singular problems is that the space decomposition must satisfy the kernel-capturing property
that is, any kernel function can be written as a sum of kernel functions drawn from the subspaces. In particular, each subspace \(V_i\) must be rich enough to support kernel functions; in our context, this is not satisfied by the choice (5.4), accounting for its poor behaviour as \(\gamma \rightarrow \infty \).
In the mesh triangulation \(\mathscr {M}_h\), we denote the star of a vertex \(v_i\) as the patch of elements sharing \(v_i\), i.e.,
This induces an associated space decomposition, called the star patch, by
This is illustrated in Fig. 1 (left). We call the induced relaxation method a star iteration. In effect, each subspace solve solves for the degrees of freedom in the interior of the patch of cells, with homogeneous Dirichlet conditions on the boundary of the patch. Given a vertex or edge midpoint \(v_i\), we denote the point-block patch \(V_i\) as the span of the basis functions associated with degrees of freedom that evaluate a function at \(v_i\) (see Fig. 1, middle). The induced relaxation method solves for all colocated degrees of freedom simultaneously. These two space decompositions coincide for the \([\mathbb {P}_1]^d\)-\(\mathbb {P}_1\) discretization.
Illustrations of the star patch of the center vertex (left) and the point-block patch (middle) for the finite element pair \([\mathbb {P}_2]^2\)-\(\mathbb {P}_1\). Note that these two patches (right) are the same for \([\mathbb {P}_1]^2\)-\(\mathbb {P}_1\) discretization. Here, black dots represent the degrees of freedom, and the blue lines gather degrees of freedom solved for simultaneously in the relaxation
We now briefly explain why these two decompositions approximately satisfy the kernel-capturing condition (5.5) for the finite element pair \([\mathbb {P}_1]^d\)-\(\mathbb {P}_1\). First, we define an approximate kernel
Since \(\mathbf {n}_k\) is the current approximation to the director \(\mathbf {n}\), we have \(\mathbf {n}_k\in V_h=\sum _i V_i\). We are therefore able to express \(\mathbf {n}_k\) as \(\mathbf {n}_k=\sum _i \mathbf {n}_k^i\), where \(\mathbf {n}_k^i\in V_i\) describes the function at the vertex \(v_i\). Similarly, we split \(\mathbf {u}_h\) into \(\mathbf {u}_h=\sum _i \mathbf {u}_h^i\) with \(\mathbf {u}_h^i\in V_i\). For each vertex \(v_i\), the requirement \(u_h \in \tilde{\mathscr {N}}_h\) yields
The definition of \(V_i\) ensures that \(\mathbf {u}_h^i\) and \(\mathbf {n}_k^i\) are only supported on the interior of the star of \(v_i\). We deduce that on each vertex
which yields \( \sum _j\mathbf {n}_k^j\cdot \mathbf {u}_h^i =\mathbf {n}_k\cdot \mathbf {u}_h^i=0\). Hence, \(\mathbf {u}_h^i\in \tilde{\mathscr {N}}_h \forall i\) and we obtain the kernel-capturing condition (5.5) for the approximate kernel \(\tilde{\mathscr {N}}_h\).
For the \([\mathbb {P}_2]^d\)-\(\mathbb {P}_1\) finite element pair, the satisfaction of the kernel-capturing property for the approximate kernel follows along similar lines. For the point-block patch, (5.7) still holds. The star patch uses larger subspaces, each one including multiple point-block patches, but it can be easily verified that (5.7) is still fulfilled.
5.1.1 Robustness analysis of the approximate kernel
While we are not able to prove the kernel capturing property for the kernel (5.2), we can still obtain the spectral inequalities
when using the approximate kernel (5.6). Here, \(D_{h,\gamma }\) is the preconditioner to be specified later for the operator \(A_{h,\gamma }\) and \(C\le D\) represents \(\Vert \mathbf {u}\Vert _C \le \Vert \mathbf {u}\Vert _D\) for all \(\mathbf {u}\). We prove that \(c_1\) depends on \(\gamma \), but the dependence can be well controlled so that the preconditioner is not badly affected by varying \(\gamma \), while \(c_2\) is always independent of \(\gamma \). For simplicity, we prove the case for the equal-constant nematic case with the \([\mathbb {P}_1]^d\)-\(\mathbb {P}_1\) discretization; extensions to the non-equal-constant cholesteric case and to the \([\mathbb {P}_2]^d\)-\(\mathbb {P}_1\) discretization are possible.
We define the operator associated to \(a^m\), \(A_{h,\gamma }: V_h\rightarrow V_h^*\), by
For the space decomposition \(V_h= \sum _i V_i\), we denote the lifting operator (the natural inclusion) by \(I_i:V_i\rightarrow V_h\) and choose the Galerkin subspace operator \(A_i: V_i\rightarrow V_i\) to satisfy
This implies that \(A_i = I_i^*A_{h,\gamma }I_i\).
The additive Schwarz preconditioner \(D_{h,\gamma }\) for a problem \(A_{h,\gamma }w_h=d_h\) associated with the space decomposition (5.3) is defined by the action of its inverse [56]:
given by
with \(w_i\in V_i\) being the unique solution of
Hence, we can rewrite the preconditioning operator \(D_{h,\gamma }^{-1}\) in operator form as
We now state for completeness a classical result in the analysis of additive Schwarz preconditioners, see e.g. [50, Theorem 3.1] and the references therein.
Theorem 5.1
Define the splitting norm for \(\mathbf {u}_h \in V_h\) as
This splitting norm is equal to the norm \(\Vert \mathbf {u}_h\Vert _{D_{h,\gamma }}:=\langle D_{h,\gamma } \mathbf {u}_h, \mathbf {u}_h\rangle _0^{1/2}\) generated by the additive Schwarz preconditioner, i.e. it holds that
To build intuition, let us examine why Jacobi relaxation defined by the space decomposition (5.4) is not robust as \(\gamma \rightarrow \infty \). With (5.4), the decomposition \(\mathbf {u}_h=\sum _i\mathbf {u}_i, \mathbf {u}_i\in V_i\) is unique. It yields that
where \(a\preceq b\) means that there exists a constant c independent of a and b such that \(a\le cb\). Note that the bound in (5.9) is parameter-dependent and deteriorates as \(\gamma \rightarrow \infty \) or \(h\rightarrow 0\).
In order to deduce the robustness result for our approximate kernel (5.6), we first derive the following lemma.
Lemma 5.1
Let \(\mathbf {u}_0=\sum _i \mathbf {u}_0^i \in \tilde{\mathscr {N}}_h\) and assume \(\mathbf {n}_k\in [\mathbb {P}_1]^d\). Then it holds that
where \(\mathrm {D}\mathbf {n}_k\) denotes the Jacobian matrix of \(\mathbf {n}_k\).
Proof
Consider the vertex \(v_i\) on the boundary of an element T. As \( \mathbf {n}_k\in [\mathbb {P}_1]^d\), we have
Note that \(\mathbf {u}_0^i \cdot \mathbf {n}_k\) vanishes at the vertex \(v_i\) as \(\mathbf {u}_0 \in \tilde{\mathscr {N}}_h\). Moreover, we know that \(\mathbf {u}_0^i(x)/\Vert \mathbf {u}_0^i(x)\Vert \) is constant on the interior of the patch around \(v_i\), and \(\mathbf {u}_0^i(x)\) is zero on the boundary of the patch, since we can write \(\mathbf {u}_0^i(x)=\mathbf {u}_0(v_i) \psi _i(x)\) with \(\psi _i\) denoting the scalar piecewise linear basis function (vanishing outside the patch) associated with \(v_i\). Therefore, we can deduce \(\mathbf {u}_0^i(x) \cdot \mathbf {n}_k(v_i)=0\) on T. In addition, we have \(\Vert x-v_i\Vert \preceq h\) on the element T. We thus conclude that
From this we are able to show that for both the star and point-block patches around \(v_i\),
Therefore, with the local support of \(\mathbf {u}_0^i\) we have
\(\square \)
We now derive the general form of the spectral bounds in (5.8). This follows a similar approach to [50, Theorem 4.1], but with a different assumption on the splitting approximation, to allow for a dependence on \(\gamma \). For brevity of notation, we respectively denote the standard \(L^2\), \(H^1\) and \(L^{\infty }\) norms by \(\Vert \cdot \Vert _0\), \(\Vert \cdot \Vert _1\) and \(\Vert \cdot \Vert _\infty \). Given a space decomposition \(V_h=\sum _i V_i\), we define its overlap \(N_O\) as
where
measures the interaction between each subspace.
Theorem 5.2
Let \(\{V_i\}\) be a subspace decomposition of \(V_h\) with overlap \(N_O\). Assume that the finite element pair \(V_h\)-\(Q_h\) is inf-sup stable for the mixed problem
where f is a known functional. Furthermore, assume that the function \(\mathbf {u}_h\in V_h\) and the kernel function \(\mathbf {u}_0\in {\mathscr {N}}_h\) can be split locally with estimates depending on the mesh size h and possibly on \(\gamma \) if the kernel-capturing property is not satisfied:
Then the additive Schwarz preconditioner \(D_{h,\gamma }\) built on the decomposition \(\{V_i\}\) satisfies
with constants \(c_1\) and \(c_2\) independent of \(\gamma \).
Proof
The upper bound can be directly given by [50, Lemma 3.2] independent of the form of partial differential equations.
For the lower bound, choose \(\mathbf {u}_h\in V_h\) and split it into \(\mathbf {u}_h = \mathbf {u}_0 + \mathbf {u}_1\), by solving
Testing with \(\mathbf {v}_h=0\) in (5.11), we obtain that
Hence, \(\mathbf {n}_k\cdot \mathbf {u}_0 = 0\) a.e., that is to say \(\mathbf {u}_0\in {\mathscr {N}}_h\).
By stability of the finite element pair \(V_h\)-\(Q_h\), we have
It implies further that
by the boundedness of \(\mathbf {n}_k\) and
by the form of the operator \(A_{h,\gamma }\), respectively. Using \(\mathbf {u}_0=\mathbf {u}_h-\mathbf {u}_1\), we have in addition that
We now calculate
completing the proof of the spectral estimates (5.10). \(\square \)
Remark 5.1
Note that in Theorem 5.2, if the kernel-capturing property (5.5) is satisfied, then \(c_3\) will be zero. Hence, we will instead get a parameter-independent result.
Corollary 5.1
In Theorem 5.2, if we take \(V_h\)-\(Q_h\) to be constructed by the \([\mathbb {P}_1]^d\)-\(\mathbb {P}_1\) element, it holds that
with constants \(c_1(h)\), \(c_2(h)\sim \mathscr {O}(h^{-2})\).
Proof
We follow the main argument of Theorem 5.2. We have only proven the kernel-capturing property for the approximate kernel (5.6) rather than (5.2), and need to account for this in the estimates. From Lemma 5.1 we have that
With the choice of \(V_h= [\mathbb {P}_1]^d\), we will use the so-called inverse inequality (its proof can be found in any finite element book, e.g., [13]) which states that
Therefore, it is straightforward to obtain that \(c_1\) and \(c_2\) are actually \(\mathscr {O}(h^{-2})\). Notice here we have also used the form of \(\Vert \cdot \Vert _{A_{h,\gamma }}\) in estimating \(c_2(h)\).
Finally, substituting the form of \(c_3\) in (5.12), we derive
with constants \(c_1(h)\), \(c_2(h)\sim \mathscr {O}(h^{-2})\). \(\square \)
The above Corollary 5.1 implies that we cannot entirely get rid of parameter \(\gamma \) in the spectral estimates if the kernel-capturing property for the modified kernel (5.2) is not satisfied and instead we get an additional factor of \(\gamma h^2 \Vert \mathrm {D}\mathbf {n}_k\Vert _\infty ^2\). However, this \(\gamma \)-dependence can be well controlled and does not impinge on the effectiveness of our smoother; the dependence improves as the mesh becomes finer or as \(\mathbf {n}_k\) becomes smoother.
5.2 Prolongation
To construct a parameter-robust multigrid method, the prolongation operator is also required to be continuous (in the energy norm associated with the PDE) with the continuity constant independent of the penalty parameter \(\gamma \) [50, Theorem 4.2]. In the context of the Oseen, Navier–Stokes, and linear elasticity equations, the prolongation operator was modified in order to guarantee that the continuity constant is \(\gamma \)-independent [11, 21, 50]. However, in our experiments with the Oseen–Frank system, we observe robust convergence with respect to \(\gamma \), even when using the (cheaper) standard prolongation (see Sect. 6.2 for specific details). This can be seen in Tables 7 and 8 of Sect. 6, for example. Hence, we will use the standard prolongation with no modification in this work.
Remark 5.2
Since both discretizations \([\mathbb {P}_1]^d\)-\(\mathbb {P}_1\) and \([\mathbb {P}_2]^d\)-\(\mathbb {P}_1\) are nested, i.e., \(V_H\subset V_h\), the standard prolongation \(P_H\) is actually a continuous (in the \(H^1\)-norm) natural inclusion.
6 Numerical experiments
6.1 Algorithm details
In the following numerical experiments, we use the \([\mathbb {P}_2]^3\)-\(\mathbb {P}_1\) element pair and use flexible GMRES [47] as the outermost linear solver, since GMRES [48] is applied in the multigrid relaxation. An absolute tolerance of \(10^{-8}\) was used for the nonlinear solver, except for the convergence rate tests in Fig. 5, which used \(10^{-10}\). A relative tolerance of \(10^{-4}\) was used for the inner linear solver. We use the full block factorization preconditioner
where \(\tilde{A}_\gamma ^{-1}\) represents solving the top-left block \(A_{\gamma }\) inexactly by our specialized multigrid algorithm and the Schur complement approximation \(\tilde{S}_\gamma ^{-1}\) is given by (3.6). The multiplier mass matrix inverse \(M_{\lambda }^{-1}\) is solved using Cholesky factorization.
For \(\tilde{A}_{\gamma }^{-1}\), we perform a multigrid V-cycle, where the problem on the coarsest grid is solved exactly by Cholesky decomposition. On each finer level, as relaxation we perform 3 GMRES iterations preconditioned by the additive star (denoted as ALMG-STAR) iteration or additive point-block Jacobi (denoted as ALMG-PBJ) iteration. In order to achieve convergence results independent of the number of cores used in parallel, we only report iteration counts using additive relaxation, although multiplicative ones generally give better convergence.
The solver described above is implemented in the Firedrake [45] library which relies on PETSc [4] for solving linear systems. The star and Vanka relaxation methods are implemented using the PCPATCH preconditioner recently included in PETSc [20].
6.2 Numerical results
All the tests are executed on a computer with an Intel(R) Xeon(R) Silver 4116 CPU@2.10GHz processor. We denote #refs and #dofs as the number of mesh refinements and degrees of freedom, respectively, in the following experiments.
6.2.1 Periodic boundary condition in a square slab
Following the nematic benchmarks in [3], we consider a generalized twist equilibrium configuration in a square \(\varOmega = [0,1]\times [0,1]\), which is proven to have an analytical solution [52]. We will investigate the robustness of the solver when applied to unequal Frank constants and nonzero cholesteric pitch.
The problem has periodic boundary conditions in the x-direction and Dirichlet boundary conditions in the y-direction, with values
where \(\theta _0={\pi }/8\).
We first consider parameter values \(K_1=1.0\), \(K_2=1.2\), \(K_3=1.0\), \(q_0=0\). The exact solution is given by
with true free energy \(2K_2\theta _0^2 \approx 0.37011\). An example of the pure twist configuration is illustrated in Fig. 2.
We use an initial guess of \(\mathbf {n}_0=[1,0,0]^\top \) in the Newton iteration and a \(10 \times 10\) mesh of triangles of negative slope as the coarse grid.
A sample solution of the twist configuration. Colors represent the magnitude of directors
We first compare in Table 1 the nonlinear convergence of the Newton linearization (2.14) against that of the Picard iteration (2.15) we propose. For these experiments we use the augmented Lagrangian preconditioner with ideal inner solvers (denoted as ALLU), i.e. where the top-left block is solved exactly by LU factorization. The Picard iteration requires substantially fewer nonlinear iterations for large \(\gamma \). We expect that this relates to the degradation of the coercivity estimate given in Lemma 2.3, and will be analyzed in future work. Similar results were obtained on other test cases and we adopt the Picard iteration henceforth.
To see the efficiency of the Schur complement approximation (3.6) we used in Sect. 3, we give the number of Krylov iterations for ALLU in Table 2. It can be observed that as \(\gamma \) increases, the preconditioner becomes a better approximation to the real Jacobian inverse and that the preconditioner is mesh-independent.
The performance of ALMG-STAR and ALMG-PBJ are illustrated in Tables 3 and 4, respectively, where both mesh-independence for \(\gamma =10^6\) and \(\gamma \)-robustness are observed.
We also test the robustness of ALMG-STAR and ALMG-PBJ on other problem parameters, the twist elastic constant \(K_2>0\) and the cholesteric pitch \(q_0\). To this end, we continue \(K_2\in [0.2, 8]\) and \(q_0\in [0,8]\) with step 0.1. We fix \(\gamma =10^6\), since it gives the best performance in Tables 3 and 4. The numerical results of ALMG-STAR and ALMG-PBJ in \(K_2\)- and \(q_0\)-continuation are shown in Figs. 3 and 4, respectively. Clearly, a stable number of linear iterations is shown for both continuation experiments.
Average number of FGMRES iterations per Newton iteration when continuing in \(K_2\) for the LC problem in a square slab
Average number of FGMRES iterations per Newton iteration when continuing in \(q_0\) for the LC problem in a square slab
To examine the convergence order of the discretization as a function of \(\gamma \), we apply the ALMG-PBJ solver for \(\gamma =10^4, 10^5\) and \(10^6\). Note that the convergence result does not rely on the solver used. Figure 5 shows the \(L^2\)- and \(H^1\)-error between the computed director and the known analytical solution. We observe third order convergence of the director in the \(L^2\) norm and second order convergence in the \(H^1\) norm for all values of \(\gamma \) considered.
The convergence of the computed director as the mesh is refined for the nematic LC problem in a square slab
To investigate the computational efficiency of the AL approach, we compare our proposed AL-based solvers (ALMG-PBJ and ALMG-STAR) with a monolithic multigrid preconditioner using Vanka relaxation [1, 53] on each level (denoted as MGVANKA) in Table 5. Essentially, MGVANKA applies multigrid to the coupled director-multiplier problem, with an additive Schwarz relaxation organised around gathering all director dofs coupled to a given multiplier dof. All results are computed in serial. In our experiments, these two AL-based solvers outperform MGVANKA even for small problems of about five thousand dofs. In particular, ALMG-PBJ is the fastest method considered and is approximately five times faster than MGVANKA for a problem with about five million dofs. We also notice that ALMG-STAR is slower than ALMG-PBJ, which is caused by the size of the star patch being larger than that of the point-block patch, requiring more work in the multigrid relaxation.
6.2.2 Equal-constant nematic case in an ellipse
Consider an ellipse of aspect ratio 3/2 with strong anchoring boundary condition \(\mathbf {n} = [0, 0, 1]^\top \) imposed on the entire boundary. We consider the equal-constant nematic case \(K_1=K_2=K_3=1\), \(q_0=0\) to verify the theoretical results presented in previous sections with corresponding discretizations. We use the initial guess \(\mathbf {n}_0 = [0, 0, 0.8]^\top \) in the nonlinear iteration. The coarsest triangulation, generated in Gmsh [26], is illustrated in Fig. 6.
The coarse mesh of the ellipse
To verify our theoretical results on the improvement of the discrete enforcement of the constraint in Sect. 4, we vary the penalty parameter \(\gamma \), use one refinement for the fine mesh, and employ the \([\mathbb {P}_1]^3\)-\(\mathbb {P}_1\) element. The data is plotted in Fig. 7. The \(L^2\)-norm \(\Vert \mathbf {n}\cdot \mathbf {n}-1\Vert _0\) of the residual of the constraint decreases as \(\gamma \) grows, and scales like \(\mathscr {O}(\gamma ^{-1/2})\) as expected.
Comparison of the computed constraint \(\Vert \mathbf {n}\cdot \mathbf {n}-1\Vert _0\) and the reference line \(\mathscr {O}(\gamma ^{-1/2})\) using the \([\mathbb {P}_1]^3\)-\(\mathbb {P}_1\) finite element pair for equal-constant nematic LC problems in an ellipse
The efficiency of the Schur complement approximation of Sect. 3 for the \([\mathbb {P}_2]^3\)-\(\mathbb {P}_1\) element can be observed in Table 6.
Tables 7 and 8 demonstrate the robustness of ALMG-STAR and ALMG-PBJ with respect to \(\gamma \) and mesh refinement for the \([\mathbb {P}_2]^3\)-\(\mathbb {P}_1\) element. It can be seen that both solvers are robust with respect to the penalty parameter \(\gamma \), and with respect to the mesh size h for \(\gamma =10^6\). The number of nonlinear iterations and the number of FGMRES iterations per Newton step remain stable.
Code availability For reproducibility, both the solver code [55] and the exact version of Firedrake used [22] to produce the numerical results of this paper have been archived on Zenodo. An installation of Firedrake with components matching those used in this paper can be obtained by following the instructions at https://www.firedrakeproject.org/download.html with
7 Conclusions
The results in this paper divide into two categories: results about the Oseen–Frank model and its discretization, and results about the augmented Lagrangian method for solving it. For the former, we extended the well-posedness results of [2] for nematic problems to the cholesteric case. We also showed that the Schur complement of the discretized system is spectrally equivalent to the Lagrange multiplier mass matrix. For the latter, we showed that the AL method improves the discrete enforcement of the constraint, and devised a parameter-robust multigrid scheme for the augmented director block. The key point in this is to capture the kernel of the semi-definite augmentation term in the multigrid relaxation. Numerical experiments validate the results and indicate that the proposed scheme outperforms existing monolithic multigrid methods.
Notes
In fact, \(\mathbf {H}_0^1(\varOmega )=\mathbf {H}_0(\mathrm {div},\varOmega )\cap \mathbf {H}_0(\mathrm {curl},\varOmega )\) holds for any bounded Lipschitz domain \(\varOmega \) [27, Lemma 2.5].
References
Adler, J.H., Atherton, T.J., Benson, T., Emerson, D.B., MacLachlan, S.P.: Energy minimization for liquid crystal equilibrium with electric and flexoelectric effects. SIAM J. Sci. Comput. 37(5), S157–S176 (2015)
Adler, J.H., Atherton, T.J., Emerson, D.B., Maclachlan, S.P.: An energy-minimization finite element approach for the Frank–Oseen model of nematic liquid crystals. SIAM J. Numer. Anal. 53(5), 2226–2254 (2015)
Adler, J.H., Atherton, T.J., Emerson, D.B., Maclachlan, S.P.: Constrained optimization for liquid crystal equilibria. SIAM J. Sci. Comput. 38(1), 50–76 (2016)
Balay, S., Abhyankar, S., Adams, M.F., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Eijkhout, V., Gropp, W.D., Kaushik, D., Knepley, M., McInnes, L.C., Rupp, K., Smith, B.F., Zhang, H.: PETSc users manual. Tech. Rep. ANL-95/11 - Revision 3.9, Argonne National Laboratory (2018)
Ball, J.M.: Mathematics and liquid crystals. Mol. Cryst. Liq. Cryst. 647(1), 1–27 (2017)
Bedford, S.J.: Calculus of variations and its application to liquid crystals. Ph.D. thesis, University of Oxford (2014)
Beik, F.P.A., Benzi, M.: Block preconditioners for saddle point systems arising from liquid crystal directors modeling. Calcolo 55(29), 1–12 (2018)
Beik, F.P.A., Benzi, M.: Iterative methods for double saddle point systems. SIAM J. Matrix Anal. Appl. 39, 902–921 (2018)
Benzi, M., Beik, F.P.A.: Uzawa-type and augmented Lagrangian methods for double saddle point systems. In: Bini, D., Di Benedetto, F., Tyrtyshnikov, E., Van Barel, M. (eds.) Structured Matrices in Numerical Linear ALgebra. Springer INdAM Series, vol. 20. Springer, Cham (2019)
Benzi, M., Golub, G.H., Liesen, J.: Numerical solution of saddle point problems. Acta Numer. 14, 1–137 (2005)
Benzi, M., Olshanskii, M.A.: An augmented Lagrangian-based approach to the Oseen problem. SIAM J. Sci. Comput. 28(6), 2095–2113 (2006)
Chandrasekhar, S.: Liq. Cryst., 2nd edn. Cambridge University Press, Cambridge (1992)
Ciarlet, P.G.: The Finite Element for Elliptic Problems. North-Holland, Amsterdam, New York, Oxford (1978)
Clément, P.: Approximation by finite element functions using local regularization. Rev. Française Automat. Informat. Recherche Opérationnelle Sér. Rouge Anal. Numér. 9(R-2), 77–84 (1975)
de Gennes, P.G.: The Physics of Liquid Crystals. Oxford University Press, Oxford (1974)
Elman, H.C., Silvester, D., Wathen, A.J.: Finite Elements and Fast Iterative Solvers: With Applications in Incompressible Fluid Dynamics, 2nd edn. Oxford University Press, Oxford (2014)
Emerson, D.B.: Advanced discretizations and multigrid methods for liquid crystal configurations. Ph.D. thesis, Tufts University (2015)
Emerson, D.B., Farrell, P.E., Adler, J.H., MacLachlan, S.P., Atherton, T.J.: Computing equilibrium states of cholesteric liquid crystals in elliptical channels with deflation algorithms. Liq. Cryst. 45(3), 341–350 (2018)
Ericksen, J.L.: Liquid crystals with variable degree of orientation. Arch. Ration. Mech. Anal. 113(2), 97–120 (1991)
Farrell, P.E., Knepley, M.G., Wechsung, F., Mitchell, L.: Pcpatch: software for the topological construction of multigrid relaxation methods. arXiv preprint arXiv:1912.08516 (2019)
Farrell, P.E., Mitchell, L., Wechsung, F.: An augmented Lagrangian preconditioner for the 3D stationary incompressible Navier-Stokes equations at high Reynolds number. SIAM J. Sci. Comput. 41, A3073–A3096 (2019)
Firedrake-Zenodo: Software used in ’Augmented Lagrangian preconditoners for the Oseen–Frank model of nematic and cholesteric liquid crystals’ (2020). https://doi.org/10.5281/zenodo.4249051
Fortin, M., Glowinski, R.: Augmented Lagrangian Methods: Applications to the Numerical Solution of Boundary-Value Problems, Studies in Mathematics and Its Applications, vol. 15. Elsevier Science Ltd, Amsterdam (1983)
Frank, F.C.: Liquid crystals. Faraday Discuss. 25, 19–28 (1958)
Friedel, V.: Les états mésomorphes de la matiére. Ann. Phys. 18, 273–474 (1922)
Geuzaine, C., Remacle, J.F.: Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. Int. J. Numer. Methods Eng. 79(11), 1309–1331 (2009)
Girault, V., Raviart, P.A.: Finite Element Methods for Navier–Stokes Equations: Theory and Algorithms, 1st edn. Springer, Cham (2011)
Glowinski, R., Le Tallec, P.: Augmented Lagrangian methods for the solution of variational problems, chap. 3. In: Studies in Applied Mathematics, pp. 45–121. SIAM (1989)
Glowinski, R., Lin, P., Pan, X.B.: An operator-splitting method for a liquid crystal model. Comput. Phys. Commun. 152(3), 242–252 (2003)
He, X., Vuik, C., Klaij, C.M.: Combining the augmented Lagrangian preconditioner with the simple Schur complement approximation. SIAM J. Sci. Comput. 40(3), A1362–A1385 (2018)
Heister, T., Rapin, G.: Efficient augmented Lagrangian-type preconditioning for the Oseen problem using Grad-Div stabilization. Int. J. Numer. Methods Fluids 71(1), 118–134 (2012)
Hu, Q., Tai, X., Winther, R.: A saddle point approach to the computation of harmonic maps. SIAM J. Numer. Anal. 47(2), 1500–1523 (2009)
John, V., Linke, A., Merdon, C., Neilan, M., Rebholz, L.: On the divergence constraint in mixed finite element methods for incompressible flows. SIAM Rev. 59(3), 492–544 (2017)
Knabner, P., Angermann, L.: Numerik partieller Differentialgleichungen. Springer, Berlin (2000)
Lee, Y., Wu, J., Xu, J., Zikatanov, L.: Robust subspace correction methods for nearly singular systems. Math. Models Methods Appl. Sci. 17(11), 1937–1963 (2007)
Lin, F.H.: Nonlinear theory of defects in nematic liquid crystals; phase transition and flow phenomena. Commun. Pure Appl. Math. 42(6), 789–814 (1989)
Lin, P., Richter, T.: An adaptive homotopy multi-grid method for molecule orientations of high dimensional liquid crystals. J. Comput. Phys. 225(2), 2069–2082 (2007)
Lin, P., Tai, X.: An Augmented Lagrangian method for the microstructure of a liquid crystal model, chap. 7. In: Fitzgibbon, W., Kuznetsov, Y., Neittaanmäki, P., Pironneau, O. (eds.) Modeling, Simulation and Optimization for Science and Technology, vol. 34, pp. 123–137. Springer, Dordrecht (2014)
Nocedal, J., Wright, S.J.: Numerical Optimization. Springer Series in Operations. Springer, Berlin (1999)
Olshanskii, M.A.: A low order Galerkin finite element method for the Navier-Stokes equations of steady incompressible flow: a stabilization issue and iterative methods. Comput. Methods Appl. Mech. Eng. 191(47), 5515–5536 (2002)
Olshanskii, M.A., Lube, G., Heister, T., Löwe, J.: Grad-div stabilization and subgrid pressure models for the incompressible Navier–Stokes equation. Comput. Methods Appl. Mech. Eng. 198(49), 3975–3988 (2009)
Oseen, C.W.: The theory of liquid crystals. Trans. Faraday Soc. 29(140), 883–899 (1933)
Polyak, V.T., Tret’yakov, N.V.: The method of penalty estimates for conditional extremum problems. USSR Comput. Math. Math. Phys. 13(1), 42–58 (1974)
Ramage, A., Gartland, E.: A preconditioned nullspace method for liquid crystal director modeling. SIAM J. Sci. Comput. 35(1), B226–B247 (2013)
Rathgeber, F., Ham, D.A., Mitchell, L., Lange, M., Luporini, F., McRae, A.T.T., Bercea, G.T., Markall, G.R., Kelly, P.H.J.: Firedrake: automating the finite element method by composing abstractions. ACM Trans. Math. Softw. 43(3), 1–27 (2017)
Reinitzer, F.: Beiträge zur Kenntnis des Cholesterins. Monatsh. Chem. 9, 421–441 (1888)
Saad, Y.: A flexible inner-outer preconditioned GMRES algorithm. SIAM J. Sci. Comput. 14(2), 461–469 (1993)
Saad, Y., Schultz, M.: GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 7(3), 856–869 (1986)
Schöberl, J.: Multigrid methods for a parameter dependent problem in primal variables. Numer. Math. 84(1), 97–119 (1999)
Schöberl, J.: Robust multigrid methods for parameter dependent problems. Ph.D. thesis, Johannes Kepler University Linz (1999)
Silvester, D., Wathen, A.J.: Fast iterative solution of stabilised Stokes systems. Part II: using general block preconditioners. SIAM J. Numer. Anal. 31, 1352–1367 (1994)
Stewart, I.W.: The Static and Dynamic Continuum Theory of Liquid Crystals: A Mathematical Introduction. CRC Press, Boca Raton (2004)
Vanka, S.P.: Block-implicit multigrid calculation of two-dimensional recirculating flows. Comput. Method. Appl. Mech. Eng. 59(1), 29–48 (1986)
Wathen, A.J., Silvester, D.: Fast iterative solution of stabilised Stokes systems. Part I: using simple diagonal preconditioners. SIAM J. Numer. Anal. 30(3), 630–649 (1991)
Xia, J.: ALpaper-numerics (2020). https://doi.org/10.5281/zenodo.4257094
Xu, J.: Iterative methods by space decomposition and subspace correction. SIAM Rev. 34(4), 581–613 (1992)
Author information
Authors and Affiliations
Corresponding author
Additional information
Communicated by Marko Huhtanen.
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This work is supported by the National University of Defense Technology, the EPSRC Centre for Doctoral Training in Partial Differential Equations [Grant Number EP/L015811/1], the EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling [Grant Number EP/L015803/1] in collaboration with London Computational Solutions, and by the Engineering and Physical Sciences Research Council [Grant Numbers EP/R029423/1 and EP/V001493/1]
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Xia, J., Farrell, P.E. & Wechsung, F. Augmented Lagrangian preconditioners for the Oseen–Frank model of nematic and cholesteric liquid crystals. Bit Numer Math 61, 607–644 (2021). https://doi.org/10.1007/s10543-020-00838-9
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10543-020-00838-9