Augmented Lagrangian preconditioners for the Oseen–Frank model of nematic and cholesteric liquid crystals

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.


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 n : Ω → 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π 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 n : Ω → 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 n and −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 n · 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-thenoptimize 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 meshindependent and the performance of the iteration can be improved by increasing the value of γ . These references also apply the discretize-then-optimize approach to tackle the pointwise unit-length vector constraint. In this paper, we will employ the optimizethen-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.

Oseen-Frank model
be an open, bounded domain with Lipschitz boundary ∂Ω and denote H 1 : v| ∂Ω = g} with a vector field g ∈ H 1/2 (∂Ω; S 2 ). Here, S 2 represents the surface of the unit ball centered at the origin. Assume that the cholesteric LC occupying the domain Ω is equipped with a rigid anchoring (Dirichlet) boundary condition n| ∂Ω = g 1 . The Oseen-Frank model (cf. [24]) considers the following minimization problem: subject to n · n = 1 a.e., (2.1) where the Frank energy density W (n) is of the form where tr(·) denotes the trace of a matrix, K i ∈ R (i = 1, 2, 3, 4) are elastic constants (called Frank constants) and q 0 ≥ 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 ∇n is matrix-valued and (∇n) 2 denotes the matrix multiplication of the matrix ∇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 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 κ = K 2 /K 3 and I is the second-order identity tensor. By the classical equality the original energy functional J (n) can be written as (2.5) Here and throughout this work, ·, · 0 denotes the inner product in L 2 (Ω) with its induced norm · 0 . It can be observed that the auxiliary tensor Z contributes to the nonlinearity of J (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 → ∞, 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 Z is uniformly (with respect to x ∈ Ω) symmetric positive definite (USPD) as long as sufficient control is maintained on n · n − 1. This property of Z plays an essential role in proving the well-posedness of the saddle-point problem in the nematic case. We restate the result of Z being USPD in the following, as it is important later:
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 4related 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 [ 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.
With a suitable spatial discretization (we only consider conforming finite elements in this work, i.e., V h ⊂ H 1 0 (Ω), Q h ⊂ L 2 (Ω)), a symmetric saddle-point system must be solved at each Newton iteration: where U and P represent the coefficient vectors of δn and δλ in terms of the basis functions of V h and Q h , respectively. We can accordingly write the discrete variational problem as: find δn h ∈ V h and δλ h ∈ Q h such that where a(·, ·) and b(·, ·) are bilinear forms given by and F and G are linear functionals in the forms of and G(μ) = − μ, n k · n k − 1 0 .

Remark 2.5
The well-posedness of the continuous and discretized Newton system (with the ([Q m ] d ⊕ B F )-Q 0 finite element pair, m ≥ 1) for a generalized nematic LC problem is discussed in [2], where B F denotes the space of quadratic bubbles and Q k represents tensor product piecewise C 0 polynomials of degree k ≥ 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 P k the set of piecewise C 0 polynomials of degree k ≥ 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 problemdependent constants. In fact, two additional q 0 -related terms in L nn 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.
Proof With the lower bound η of Z, we compute the bilinear form: where the first inequality comes from the assumption that λ 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 where ν is the outward unit normal on the boundary ∂Ω. Then using the classical It follows that 0 (the positivity follows from the assumptions) and α 0 = c/C 2 , we find that the coercivity (2.12) holds. In particular, when κ = 1 (i.e., K 2 = K 3 ), we have Z = I and thus η = 1. Then, the bilinear form becomes .
So far, the coercivity of the bilinear form a(·, ·) has been shown for all functions in H 1 0 (Ω). The discrete coercivity follows if a conforming finite element for the director space is chosen.
The boundedness of the bilinear form a(·, ·) and the right hand side functionals F(·) and G(·) 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(·, ·) for a finite element pair V h -Q h , i.e. whether there exists a constant c such that  [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 [P 2 ] 2 -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(·, ·) 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Ã −1 andS −1 for the director block and the Schur complement S = −B A −1 B , 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.

Augmented Lagrangian form
Now, we employ the AL stabilization strategy and modify the linearized saddle point system to control its Schur complement S.

Penalizing the constraint
We penalize the continuous form of the nonlinear constraint n · n = 1 in the AL algorithm and obtain the Lagrangiañ for γ ≥ 0. The weak form of the associated first-order optimality conditions is to find The Newton linearization at a given approximation (n k , λ k ) yields a system of the form: Thus, we have to solve the augmented discrete variational problem: where Comparing (2.14) to the original system (2.11), only the bilinear form a(·, ·) and the right hand side functional F(·) have changed. The boundedness of F c (·) follows straightforwardly via the Cauchy-Schwarz inequality. As for the coercivity of a c (·, ·), an additional assumption on the penalty parameter γ is needed.
The condition α 0 > 2γ |α − 1| in Lemma 2.3 indicates a limit on the value of γ to ensure the solvability of the augmented system (2.14). However, it is desirable to use large values of γ 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γ n k ·n k −1, v·v 0 from the linearized equations. This yields the linearized problem: 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 γ . 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 n k · u, n k · v 0 and l denotes the assembly of −2 n k · n k − 1, n k · v 0 . Note that compared to the non-augmented version (2.10), the (1, 1) block in (2.17) has an additional semi-definite term γ A * with a large coefficient γ . 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.

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.
We define the fluctuation operator κ := I − Π Q h where I : Q → Q is the identity mapping. Therefore, one has For u h , v h ∈ V h , one can split the term 4 n k · u h , n k · v 0 into the following terms using the properties of κ and Π Q h : Note here that the assembly of the first term is B M −1 λ B, where M λ 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 κ(2n k · u), κ(2n k · v) 0 characterizes the difference between A * and B M −1 λ B. The next result shows that it vanishes as the mesh size h → 0 (see Theorem 3.1) and thus, in this limit, the tractable term B M −1 λ B dominates A * . Theorem 3.1 Let (δn h , δλ h ) ∈ V h × Q h be the solution of the augmented discrete system (2.15) with corresponding degrees of freedom (U , P) ∈ R n × R m . Then, for the Newton linearization at a given approximation (n k , λ k ) satisfying α ≤ |n k | 2 ≤ β with 0 < α ≤ 1 ≤ β and |∇n k | bounded pointwise a.e., we have where · R n denotes the Euclidean norm.

Proof
Assuming v h ∈ V h and using the basis representations in V h = span{φ i } for δn h and v h : by applying the Cauchy-Schwarz inequality. One readily sees that G 1 ≤ C 1 for a certain constant C 1 from the continuity of κ. 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 v h : for some C r > 0. Then with the fact that 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 n k , ∇n k a.e. and the fact that Combining these estimates regarding G 1 , G 2 , G 3 , we find 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 −1 γ , i.e., solving linear systems involving S γ . For large γ , 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 −γ −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 −γ −1 M λ will be an effective approximation for S γ , with the approximation improving as γ → ∞.
In fact, the approximation of the inverse of the discretely augmented Schur complement (3.2) can be improved further by combining −γ M −1 λ with a good approximation of the unaugmented Schur complement S [30]. Given an approximationS 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 B A −1 B is spectrally equivalent to M λ for the equal-constant nematic case. This gives us good insight into the choice ofS −1 . Proof For the equal-constant model with Dirichlet boundary conditions n = g ∈ H 1/2 (∂Ω; S 2 ), its corresponding Lagrangian is After Newton linearization and introducing conforming finite dimensional spaces where n k and λ k represent the current approximations to n and λ, respectively. This can be rewritten in block matrix form as where as before U ∈ R n and P ∈ 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 ∇δn h , ∇v h 0 + 2 λ k , δn h · v h 0 . The coercivity property of the bilinear form from Lemma 2.2 ensures that A is positive definite. The coefficient matrix A is symmetric and indefinite (resulting in A possessing both positive and negative eigenvalues). Moreover, 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 nonnegativity of λ 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 Thus, we have where the maximum is attained at z = (P B A −1/2 ) . 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 B A −1/2 ) . This gives rise to Therefore for inf-sup stable finite element pairs, we have by (3.4) and (3.5) This indicates that B A −1 B is spectrally equivalent to M λ .

Remark 3.1
It follows from Theorem 3.2 that γ = 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 γ = 0). One should also notice that such mesh-independence for γ = 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 −1 γ is given by

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 γ . 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 γ 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 n| ∂Ω = g, where g is a nonzero constant vector satisfying |g| = 1. We use the [P 1 ] d -P 1 finite element pair in this section, so both the director n and the Lagrange multiplier λ 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 ∩ H 1 g (Ω) and Q h ⊂ L 2 (Ω), respectively, and denote V h,0 = V h ∩ H 1 0 (Ω). We restate the associated nonlinear discrete variational problem as follows: find (4.2) Note that in this step we have used the fact that since g is a constant vector, its derivative is zero.
As (4.1b) is valid for arbitrary μ h ∈ Q h and one can easily verify that n h · g ∈ Q h , we have n h · g, n h · n h − 1 0 = 0. Then taking μ h = 1 and μ h = λ h leads to 1, n h · n h − 1 0 = 0 and λ h , n h · n h − 1 0 = 0, 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 |Ω| denotes the measure of the domain Ω.
Hence, by combining (4.4) and (4.5), we have Since the right hand side of (4.6) is a fixed constant independent of γ , taking γ larger value forces the constraint approximation error n h · n h − 1 0 to become smaller. In fact, (4.6) implies that n h · n h − 1 0 ≤ O(γ −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.

A robust multigrid method for A
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 γ , as it adds a semi-definite term with a large coefficient. For the augmented Lagrangian strategy to be successful, we require a γ -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 γ . 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 equalconstant 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 ,γ : V H → V * H corresponds to the partial differential equations (PDEs) on V H . For the domain Ω, we consider a non-overlapping triangulation M H , i.e., The fine grid M h with h = H /2 is obtained by a regular refinement of the simplices in M H . In what follows we consider both the [P 1 ] d -P 1 and [P 2 ] d -P 1 discretizations.

Relaxation
After applying the AL method introduced in Sect. 2.2.1, the discrete linear variational form corresponding to the top-left block A γ = A + γ A * is given by being the trial function and v h ∈ V h the test function. Note that n k and λ k are the current approximations to the director n and the Lagrange multiplier λ, respectively, in the Newton iteration. The first two terms of a m are symmetric and coercive because of the uniform non-negativity of λ k in the assumption of our well-posedness result. The kernel of the semi-definite term involving γ is In the case of γ 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 = span(φ 1 , . . . , φ N ), Jacobi and Gauss-Seidel are induced by the space decomposition V i = span(φ i ), (5.4) 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 parameterrobustness 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 γ → ∞.
In the mesh triangulation 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 [P 1 ] d -P 1 discretization.
We now briefly explain why these two decompositions approximately satisfy the kernel-capturing condition (5.5) for the finite element pair [P 1 ] d -P 1 . First, we define an approximate kernel (5.6) Fig. 1 Illustrations of the star patch of the center vertex (left) and the point-block patch (middle) for the finite element pair [P 2 ] 2 -P 1 . Note that these two patches (right) are the same for [P 1 ] 2 -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 Since n k is the current approximation to the director n, we have n k ∈ V h = i V i . We are therefore able to express n k as The definition of V i ensures that u i h and n i k are only supported on the interior of the star of v i . We deduce that on each vertex which yields j n j k · u i h = n k · u i h = 0. Hence, u i h ∈Ñ h ∀i and we obtain the kernel-capturing condition (5.5) for the approximate kernelÑ h .
For the [P 2 ] d -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 pointblock patches, but it can be easily verified that (5.7) is still fulfilled.

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,γ is the preconditioner to be specified later for the operator A h,γ and C ≤ D represents u C ≤ u D for all u. We prove that c 1 depends on γ , but the dependence can be well controlled so that the preconditioner is not badly affected by varying γ , while c 2 is always independent of γ . For simplicity, we prove the case for the equal-constant nematic case with the [P 1 ] d -P 1 discretization; extensions to the non-equal-constant cholesteric case and to the [P 2 ] d -P 1 discretization are possible. We define the operator associated to a m , A h,γ : For the space decomposition V h = i V i , we denote the lifting operator (the natural inclusion) by I i : V i → V h and choose the Galerkin subspace operator A i : This implies that A i = I * i A h,γ I i . The additive Schwarz preconditioner D h,γ for a problem A h,γ w h = d h associated with the space decomposition (5.3) is defined by the action of its inverse [56]: Hence, we can rewrite the preconditioning operator D −1 h,γ 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 u h ∈ V h as
This splitting norm is equal to the norm u h D h,γ := D h,γ u h , u h 1/2 0 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 γ → ∞. With (5.4), the decomposition u h = i u i , u i ∈ V i is unique. It yields that where a b means that there exists a constant c independent of a and b such that a ≤ cb. Note that the bound in (5.9) is parameter-dependent and deteriorates as γ → ∞ or h → 0.
In order to deduce the robustness result for our approximate kernel (5.6), we first derive the following lemma.
where Dn k denotes the Jacobian matrix of n k .
Proof Consider the vertex v i on the boundary of an element T . As n k ∈ [P 1 ] d , we have Note that u i 0 · n k vanishes at the vertex v i as u 0 ∈Ñ h . Moreover, we know that u i 0 (x)/ u i 0 (x) is constant on the interior of the patch around v i , and u i 0 (x) is zero on the boundary of the patch, since we can write u i 0 (x) = u 0 (v i )ψ i (x) with ψ i denoting the scalar piecewise linear basis function (vanishing outside the patch) associated with v i . Therefore, we can deduce u i 0 (x) · n k (v i ) = 0 on T . In addition, we have x − v i 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 u i 0 we have 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 γ . For brevity of notation, we respectively denote the standard L 2 , H 1 and L ∞ norms by · 0 , · 1 and · ∞ . Given a space decomposition V h = 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 u h ∈ V h and the kernel function u 0 ∈ N h can be split locally with estimates depending on the mesh size h and possibly on γ if the kernel-capturing property is not satisfied:

Then the additive Schwarz preconditioner D h,γ built on the decomposition {V
with constants c 1 and c 2 independent of γ .
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 u h ∈ V h and split it into u h = u 0 + u 1 , by solving Testing with v h = 0 in (5.11), we obtain that Hence, n k · u 0 = 0 a.e., that is to say u 0 ∈ N h . By stability of the finite element pair V h -Q h , we have It implies further that u 1 1 u h 0 by the boundedness of n k and by the form of the operator A h,γ , respectively. Using u 0 = u h −u 1 , we have in addition that u 0 1 u h 1 .
We now calculate completing the proof of the spectral estimates (5.10).

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.
Proof We follow the main argument of Theorem 5.2. We have only proven the kernelcapturing 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 = [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 O(h −2 ). Notice here we have also used the form of · A h,γ in estimating c 2 (h).
Finally, substituting the form of c 3 in (5.12), we derive The above Corollary 5.1 implies that we cannot entirely get rid of parameter γ 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 γ h 2 Dn k 2 ∞ . However, this γ -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 n k becomes smoother.

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 γ [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 γ -independent [11,21,50]. However, in our experiments with the Oseen-Frank system, we observe robust convergence with respect to γ , 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 [P 1 ] d -P 1 and [P 2 ] d -P 1 are nested, i.e., V H ⊂ V h , the standard prolongation P H is actually a continuous (in the H 1 -norm) natural inclusion.

Algorithm details
In the following numerical experiments, we use the [P 2 ] 3 -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Ã −1 γ represents solving the top-left block A γ inexactly by our specialized multigrid algorithm and the Schur complement approximationS −1 γ is given by (3.6). The multiplier mass matrix inverse M −1 λ is solved using Cholesky factorization. ForÃ −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].

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.

Periodic boundary condition in a square slab
Following the nematic benchmarks in [3], we consider a generalized twist equilibrium configuration in a square Ω = [0, 1] × [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.
We use an initial guess of n 0 = [1, 0, 0] in the Newton iteration and a 10 × 10 mesh of triangles of negative slope as the coarse grid.
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 γ . 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.
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,    To examine the convergence order of the discretization as a function of γ , we apply the ALMG-PBJ solver for γ = 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 γ considered.
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 The coarse mesh of the ellipse

Equal-constant nematic case in an ellipse
Consider an ellipse of aspect ratio 3/2 with strong anchoring boundary condition n = [0, 0, 1] 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 n 0 = [0, 0, 0.8] in the nonlinear iteration. The coarsest triangulation, generated in Gmsh [26], is illustrated in Fig. 6.
To verify our theoretical results on the improvement of the discrete enforcement of the constraint in Sect. 4, we vary the penalty parameter γ , use one refinement for the fine mesh, and employ the [P 1 ] 3 -P 1 element. The data is plotted in Fig. 7. The L 2 -norm n · n − 1 0 of the residual of the constraint decreases as γ grows, and scales like O(γ −1/2 ) as expected.
The efficiency of the Schur complement approximation of Sect. 3 for the [P 2 ] 3 -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 γ and mesh refinement for the [P 2 ] 3 -P 1 element. It can be seen that both solvers are robust with respect to the penalty parameter γ , and with respect to the mesh size h for γ = 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     (5) 11.60 (5) 6.80 (5) 2.50 (6) paper can be obtained by following the instructions at https://www.firedrakeproject. org/download.html with python3 firedrake-install---doi 10.5281/zenodo.4249051

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.
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/.