High order numerical simulations for the polymer self-consistent field theory using the adaptive virtual element and spectral deferred correction methods

by   Kai Jiang, et al.

In this paper, we propose an efficient and accurate numerical method, which combine the high order adaptive virtual element method (VEM) and spectral deferred correction (SDC) scheme together, to simulate the self-consistent field theory (SCFT) model on arbitrary domains. The VEM is very flexible in handling rather general polygonal elements, which is more suitable to discrete the model defined on complex geometric domain, and the SDC scheme can efficiently improve the approximation order of contour derivative. Moreover, an adaptive method equipped with a new marking strategy is developed to efficiently solve the SCFT for strong segregation systems. Results indicate that our approaches can efficiently produce high-accuracy numerical results.


page 10

page 13

page 14

page 15

page 17


An adaptive high-order surface finite element method for the self-consistent field theory on general curved surfaces

In this paper, we develop an adaptive high-order surface finite element ...

Weakly imposed Dirichlet boundary conditions for 2D and 3D Virtual Elements

In the framework of virtual element discretizazions, we address the prob...

High-order energy stable schemes of incommensurate phase-field crystal model

This article focuses on the development of high-order energy stable sche...

High order interpolatory Serendipity Virtual Element Method for semilinear parabolic problems

We propose an efficient method for the numerical solution of a general c...

A lowest-order locking-free nonconforming virtual element method based on the reduced integration technique for linear elasticity problems

We develop a lowest-order nonconforming virtual element method for plana...

Adaptive Hermite Spectral Methods in Unbounded Domains

Recently, new adaptive techniques were developed that greatly improved t...

1 Introduction

Block polymers have attracted tremendous attention for many years due to their industrial applications that rely on the customized microstructures. In practical environment, geometric restriction strongly influences the formation of microstructures, which also provides a new opportunity to engineer novel structures. Concretely speaking, the confining geometries, resulting in structural frustration, confinement-induced entropy loss, and surface interactions, can lead to novel morphologies which are not obtained in the bulk. There are many industrial applications for the block-copolymer ordered structures at the nanoscale, such as the construction of high-capacity data storage devices, waveguides, quantum dot arrays, dielectric mirrors, nanoporous membranes, nanowires, and interference lithography charlotte2011interplay ; segalman2005patterning .

Modeling and numerical simulation provide an effective means to investigate the phase separation behavior of polymer systems. Fully atomistic and coarse-graining approaches are both computational intensive methods for calculating equilibrium microstructures of polymer systems, especially for larger and more complicated geometries wang2015molecular ; sethuraman2014 . A more and effective continuum approach is the self-consistent field theory (SCFT), which is one of the most successful modern tools for studying the phase behaviors of inhomogeneous polymer systems, such as self-assembly, thermodynamic stability. SCFT can efficiently describe polymer architecture, molecular composition, polydispersity, polymer subchain types, interaction potential and related information as a series of parameters. SCFT modeling is started with a coarse grained chain and microscopic interaction potentials used in particle model, then transforms the particle-based model into a field-theoretic framework, finally obtains a mean-field equations system within saddle-point approximation fredrickson2006equilibrium .

From the viewpoint of mathematics, the SCFT model is a complicated variational problem with many challenges, such as saddle-point, nonlinearity, multi-solutions, and multi-parameters. It is difficult to obtain an analytical solution for this model. Numerical simulation is a feasible tool to study of the polymer SCFT, which usually contains four parts: screening initial values xu2013strategy ; jiang2010spectral ; jiang2013discovery

, solving time-dependent partial differential equations (PDEs) 

ouaknin2016scft ; ceniceros2019efficient , evaluating (monomer) density operators ceniceros2019efficient , and finding saddle-points ceniceros2004numerical ; thompson2004improved ; jiang2015analytic . The equilibrium state solution of the SCFT corresponds to an ordered microstructure. Due to the delicate energy difference among different ordered patterns in polymer systems, a high order numerical method is strongly needed.

In the past several decades, spectral methods, especially the Fourier spectral method, have been the predominant tool for solving SCFT model matsen1994stable ; rasmussen2002improved ; cochran2006stability . This approach has high order precision and is efficient if a spectral collocation method can be found. However, the spectral method uses the global basis functions to discrete the spatial functions, which limits its applications on the model defined on complex geometric domains and complex boundary conditions. An alternative approach is using local basis functions to discretize spatial functions, such as the finite element method (FEM) ackerman2017finite ; wei2019finite . The precision of the FEM depends on the size and quality of the mesh and the order of local polynomial basis functions. Combined with adaptive method binev2004 , FEM can obtain better numerical accuracy with less calculation cost. However, there exist some inconvenience in the use of the FEM with adaptive method, especially when the adaptive mesh contains hanging nodes carstensen2009 , polygonal or concave mesh elements.

To address these problems, in this work, we develop a novel approach to solve SCFT model based on the virtual element method (VEM) B2013vem ; L2013vem ; lfal2013vem and spectral deferred correction (SDC) method. The VEM can be considered as an extension of conforming finite element methods to polyonal meshes, which has been developed for solving a variety of partial differential equations, see Antonietti2014vem ; Zhao2016vem ; Chen2017vem and references therein. The SDC method is a high order time discrete scheme developed recently, see ceniceros2019efficient ; Dutt2000sdc . Our contribution in this paper contains: (a) formulating the SCFT problem in real space using a high order VEM based variational form, (b) incorporating the SDC scheme into VEM for polymer systems which guarantees the high-order precision for the contour variable, (c) proposing a new adaptive VEM to solve the SCFT model, especially for highly segregated systems.

The remaining sections are organized as follows. In Sec. 2 we give the SCFT model defined on the general domain using the diblock Gaussian chains as an example. In Sec. 3, we present our numerical approach, including the (adaptive) VEM and the SDC method to solve SCFT. In Sec. 4, we demonstrate the precision and the efficiency of our method by several numerical experiments. In the Sec. 5, we summarize this work.

2 Self-consistent field theory

In this section, we give a brief introduction to the SCFT model for an incompressible AB diblock copolymer melt on the general domain . We consider a system with conformationally symmetric diblock copolymers and each has and arms joined together with a covalent bond. The total degree of polymerization of a diblock copolymer is , the -monomer fraction is , and the -monomer fraction is . The field-based Hamiltonian within mean-field approximation for the incompressible diblock copolymer melt is fredrickson2006equilibrium ; cochran2006stability


where is the Flory-Huggins parameter to describe the interaction between segments and . The terms and can be viewed as fluctuating pressure and exchange chemical potential fields, respectively. The pressure field enforces the local incompressibility, while the exchange chemical potential is conjugate to the difference of density operators. is the single chain partition function, which can be computed according to


The forward propagator

represents the probability weight that the chain of contour length

has its end at position . The variable is used to parameterize each copolymer chain such that represents the tail of the block and is the junction between the and blocks. According the flexible Gaussian chain model fredrickson2006equilibrium , satisfies the following PDE


with the initial condition and being the radius of gyration. The above PDE is well-defined through possessing an appropriate boundary condition. In this work, we consider the homogeneous Neumann boundary condition


The reverse propagator , which represents the probability weight from to , satisfies Eqn. (3) only with the right-hand side of Eqn.(3a) multiplied by . The initial condition is . The normalized segment density operators and follow from functional derivatives of with respect to and and the familiar factorization property of propagators


The first variations of the Hamiltonian with respect to fields and lead to the mean-field equations


The equilibrium state, i.e., , of the SCFT model corresponds to the ordered structure. Within the standard framework of SCFT, finding the stationary states requires the self-consistent iterative procedure as shown in the following flowchart.

Given an arbitrary domain and initial fields ,

Calculate propagators and

Compute , density operators and , and evaluate the Hamiltonian

Update fields and

Is Hamilton difference less than a prescribed tolerance ?

Converged result



Flowchart of SCFT iteration.

The propagator equation is dependent on the potential fields, and . In order to start the process, the values of and must be initialized. If the initial values are homogeneous, the gradient term in the modified diffusion equation goes to zero, leaving no driving force for the formation of a microstructure. To prevent this, there must be some spatial inhomogeneity in the initial values. For a targeted periodic structure, using the space group symmetry is a useful strategy to screen the initial configuration jiang2010spectral ; jiang2013discovery . Once initial values are ready, high-accuracy numerical methods to solve the propagator equation, and evaluate the density functions, are required to solve the SCFT model, which is also the main work in this paper. We will detail our approach in the Sec. 3.

The iteration method to update fields is dependent on the mathematical structure of SCFT. An important fact is that the effective Hamiltonian (1) of diblock copolymers can reach its local minima along the exchange chemical field , and achieve the maxima along the pressure field fredrickson2006equilibrium . Thus alternative direction gradient approaches, such as the explicit Euler method, can be used to find the saddle point. In particular, the explicit Euler approach is expressed as


An accelerated semi-implicit scheme has been developed to find the equilibrium states ceniceros2004numerical ; jiang2015analytic

. However, the existing semi-implicit method is based on the asymptotic expansion and global Fourier transformation and can not be straightforwardly applied to the local basis discretization schemes.

3 Numerical methods

Solving the propagator equations is the most time-consuming part of the entire numerical simulation, and we will discuss its implementation in detail in this section, which include the (adaptive) VEM and the SDC method. In the following, we use the to denote the common norm over a finite domain .

3.1 VEM discretization for the spatial variable

VEM is a generalization of the finite element method that is inspired from the modern mimetic finite difference scheme L2013vem . Compared with FEM, VEM can handle very general (even non-convex) polygonal elements with an arbitrary number of edges. Furthermore, VEM can naturally treat the handing nodes appearing in the mesh adaptive process as the vertices of the polygonal elements, which greatly simplifies the design and implementation of mesh adaptive algorithms. Fig. 1 gives a schematic mesh which the VEM can deal with.

Figure 1: A schematic mesh of the VEM including hanging nodes and concave polygons.

Subsequently, we will introduce the virtual element space, and discretize propagator equations (3) based on the variational formulation.

3.1.1 Virtual element space

Let be the polygonal decomposition of a given domain including a finite number of non-overlapping polygons. For any polygon element , let be the set of all edges , edges of each element, the centroid, the diameter, and the area of the element . Let be the polynomials space of degree up to on , , and be the scaled monomial basis set of with form lfal2013vem


and . We will also use instead of , where is a one-dimensional index of the nature correspondence of , for example,


The local virtual element space can be defined as B2013vem ; lfal2013vem


where denotes the common Laplace operator. is a set of polynomials of degree up to on . The dimension of is


where is the number of vertices of . The function can be defined through satisfying the following three conditions:

  • is a polynomial of degree on each edge ;

  • is globally continuous on ;

  • is a polynomial of degree in .

Correspondingly, the degree of the freedom of the contains:

  • the value of at the vertices of ;

  • the value of at the internal Gauss-Lobatto quadrature points on e;

  • the moments up to order

    of in : .

Then the global virtual element space can be defined based on the local space


The dimension of is


where , and are the total number of vertices, edges, and elements of , respectively. Since is a separable Hilbert space, it can give a set of basis functions for such that, for each


where is coefficient of the degree of the freedom corresponding to . It should be emphasized that the basis functions in the VEM do not have explicit expression as the FEM has. In practical implementation, the quantities related to the basis functions can be obtained through polynomial space, as discussed in the subsequent sections. We name order VEM by the order of polynomial space . For example, the linear and quadrature VEMs mean that the order of the projection polynomial space is and , respectively.

3.1.2 Variational formulation

Using VEM to solve PDEs (3) is based on the variational formulation whose continuous version is: find such that, for all ,


where represents the inner product. In numerical computation, the spatial function must be discretized in the finite dimensional virtual element space . Then the continuous variational formulation (17) is discretized as: find such that


Let , using the expression (16), . The discretized variational formulation (18) has the matrix form





The stiffness matrix , the mass matrix , and the cross mass matrix can be obtained through projecting local virtual element space onto polynomial space. In the sequential subsections, we will present the construction methods for local stiffness, mass, and cross mass matrices. The corresponding global matrices , and can be obtained as the standard assembly process of FEM once we have the local ones.

3.1.3 Construction of the stiffness matrix

The stiffness matrix in the VEM can be computed by the local project operator ,


which projects the local virtual element space onto the polynomial space with degree up to . For each , we have the orthogonality condition


The above condition defines only up to a constant. It can be fixed by prescribing a projection operator onto constants requiring


can be chosen as


where is the number of vertices of .

Next we compute local stiffness matrix on the polygon ,


With the projector , can be splitted into


then Eqn. (25) becomes


Replacing the second term as

we can obtain the approximate local stiffness matrix


3.1.4 Construction of the mass matrix

The mass matrix in the VEM can be obtained from the local projection, . For each ,


can not be calculated directly. Next, we show how to compute the local mass matrix L2013vem


Similar to the construction method of stiffness matrix, we can define the basis function through projection operator




Replacing the second term in the above equation as

the local mass matrix can be approximated as


3.1.5 Cross mass matrix

The local cross mass matrix on can be defined as


Applying the projection , as defined in the above Sec. 3.1.4, into the cross term, the local mass matrix can be calculated as


3.1.6 Space integral

In this section, we present the integration approach over an arbitrary polygon . We divide the polygon into triangles by linking two endpoints of each edge and the barycenter. Then we apply the common Gaussian quadrature in each triangle, and summarize these integration values.


where is the set of quadrature points of , and the corresponding quadrature weights.

3.2 SDC method for the contour variable

The deferred correction approach Dutt2000sdc first solves the PDE with an appropriate method, then uses the residual equation to improve the approximation order of numerical solution. The key idea of SDC is to use a spectral quadrature ceniceros2019efficient , such as a Gaussian or a Chebyshev-node interpolatory quadrature, to integrate the contour derivative, which can achieve a high-accuracy numerical solution with a largely reduced number of quadrature points. The detail will be presented in the sequential content.

In this work, we use the variable step Crank-Nicholson (CN) scheme to solve the semi-discretize propagator equation (19), and obtain the initial numerical solution .


where is the time step size, () is the Chebyshev node clen1960integral . It should be pointed out that other stable time schemes can be employed to solve semi-discretize propagator equation (19), such as second-order operator-splitting method rasmussen2002improved , implicit-explicit Runge-Kutta scheme ascher1997rungekutta .

Then we use the deferred correction scheme to achieve a high-accuracy numerical solution. We can give the exact semi-discretize solution of propagator by integrating (19) along the contour variable


The error between the numerical solution and exact semi-discretize solution is defined as


Multiplying both sides by , we have


the residual


can be compute by the spectral integral method with Chebyshev-nodes as presented in the Sec. 3.3. By the definition of residual , we have the error integration equation


Taking the first derivative of the above equation with respect to leads to


which can be also solved by the CN scheme (37). Then the corrected numerical solution is


Repeating the above process, one can have , is the pre-determined number of deferred corrections. The convergent order of deferred correction solution along the contour parameter is


where , is the order of the chosen numerical scheme to solve Eqns. (19) and (46). For the CN scheme, .

In summary, under appropriate regularity hypothesis, one can prove the estimator for the numerical solution



is the true solution of propagator, and .

3.3 Spectral integral method along the contour variable

In this section, we discuss the Chebyshev-node interpolatory quadrature method to integrate the residual error of Eqn. (44) for the contour variable , which has spectral accuracy for smooth integrand trefethen2008guass . The proposed scheme can be also applied to the evaluate the density operators (5) and (6). These problems can be summarized to the following integral


the integrand is a smooth function. After changing variables, the general integral becomes



. The interpolate polynomial of

at Chebyshev nodes , .




. In practice, coefficients are calculated by the fast discrete cosine transform.


Due to , we have


3.4 Adaptive VEM

The adaptive method is an important technique to improve accuracy of the solution and reduce the computational complexity. The following is the adaptive process used in SCFT calculation:

Step 1

Solve the SCFT model and obtain the numerical solution on the current mesh.

Step 2

Estimate error on each element from current numerical results.

Step 3

Mark mesh elements according to the error estimate.

Step 4

Refine or coarsen the marked elements.

Next we present some implementation details of the above adatpive process.

The estimator is an important part of the adaptive method. Let be the error of indicator function over each element ,


is the harmonic average operator huang2012


is the number of elements with as a vertex. The indicator function is an essential part in adaptive methods. In SCFT model, several spatial functions can be used as indicator functions, such as field functions, density functions, and propagators. To choose efficient indicator function, we observe the distribution of these spatial functions when the SCFT calculation converges. As an example, Fig. 2 presents the equilibrium states of , and propagator function of the last contour point , respectively, with . As one can see, the distributions of three spatial functions are similar, however, has the sharpest interface. If the numerical error of can be reduced through the adaptive method, the error of other spatial functions obviously reduce with it. Therefore, in the current adaptive method, we choose as the indicator function in the posterior error estimator.

Figure 2: The equilibrium distributions of , , and when .

Given an effective and reliable posterior error estimator , it is required a marking strategy to mark mesh elements. Classical marking strategies such as the maximum jarausch1986 and the criterion dorfler1996 , usually refine or coarsen marked mesh elements one time in one adaptive process. It may make less use of the information of posteriori error estimator. To improve it, we propose a new marking strategy, named criterion, as following


where is a positive constant, is the mean value of all element estimator , and is the nearest integer function. , and represent that cell is unchanged, refined times, and coarsened times, respectively. Obviously, this new marking criterion not only denotes which mesh element needs to be improved, but also provides the times of refinement or coarseness. Due to the limitation of the topic in this paper, we will discuss the marking criterion in detail in our further paper wei2020new .

In practical implementation, we choose a widely used quadtree approach to refine and coarsen the mesh in our SCFT simulation. Quadtree is a hierarchical data structure commonly used in adaptive numerical simulation osher2006level . Given an initial (maybe unstructured) quadrilateral mesh, quadtree approach takes its every quadrilateral element as a root node of a quadtree, then refines and coarsens every quadtree according to the above posteriori error estimator and marking criterion. Fig. 3 shows the refinement and coarsen process based on the qaudtree structure. Notice that hanging nodes in quadrilateral meshes (e.g. node in Fig. 3) in VEM are treated as polygon vertices and do not be specially treated. Therefore adaptive VEM is implemented on general polygonal meshes.

Figure 3: The adaptive mesh in VEM based on quadtree structure. (a) Refine the root 0 into 4 child nodes. (b) Refine nodes 1 and 3 into 4 child nodes, respectively. (c) Coarsen child nodes 9, 10, 11, and 12 to node 3.

4 Numerical results

In the following numerical examples, we use linear () and quadratic () VEMs to discretize the spatial variable. Due to the limitation of spatial discretization order, in time direction, we just correct the initial numerical solution one time in the SDC scheme. All the numerical examples are implemented based on the FEALPy package fealpy2020 .

4.1 The efficiency of the proposed numerical method

In this subsection, we demonstrate the effectiveness of the proposed numerical method by solving a parabolic equation and SCFT model.

4.1.1 The efficiency of solving a parabolic equation

We apply our apporach to the following parabolic equation (59)


with exact solution .

First we verify the convergent order of the linear and quadratic VEMs. The time discretization is the CN scheme using to guarantee the enough time discretization accuracy. Tab. 1 gives the error and convergent order of VEM which is consistent with the theoretical result.

Nodes Linear VEM Quadratic VEM
order order
289 4.7737e-02 1.2062e-03
1089 1.3267e-02 1.84 1.5084e-04 2.99
4225 3.4013e-03 1.96 1.8863e-05 2.99
16641 8.5563e-04 1.99 2.3582e-06 3.00
Table 1: The error order of VEM.

Second we verify the error order of the CN and SDC schemes for solving (59). For the SDC scheme, we obtain a new solution by correcting the initial numerical solution calculated by the CN scheme just once. For the spatial discretization, we use the quadratic VEM with nodes to guarantee the spatial discretization accuracy. Tab. 2 gives the convergent order of the time discretization schemes which are also consistent with theoretical results. Notice that the error showed above is the error between the true solution and the VEM numerical solugion at the end time .

order order
4 6.0605e-03 5.7514e-04
8 1.5074e-03 2.00 1.0163e-05 5.82
16 3.7637e-04 2.00 6.4626e-07 3.97
32 9.4065e-05 2.00 4.2283e-08 3.94
Table 2: The error order of the time discretization schemes.

Third we verify the integral accuracy of numerical solution along the contour variable which is required in many places of SCFT simulation, including solving PDEs and evaluating density functions. We use the quadratic VEM (66049 nodes) to discretize the parabolic equation (59) to obtain a semi-discrete system. Correspondingly the exact solution of (59) can be discretized into . Then we solve the semi-discretize system using the CN and the SDC schemes for to obtain the numerical solutions and , respectively. We integrate and along from to using a modified fourth order integral scheme press1992numerical and the spectral integral quadrature as discussed in the Sec. 3.3, respectively. The integrated values are denoted by and . The exact integral about along from to can be obtained as . The error is defined as


where . Tab. 3 presents the numerical results, and one can find that can achieve about only requiring contour discretized nodes, while requiring nodes.

4 4.2343e-03 2.4219e-04
8 1.0728e-03 4.4032e-06
16 2.6970e-04 4.0527e-06
32 6.7675e-05 4.0519e-06
64 1.7401e-05 4.0520e-06
128 5.8778e-06 4.0520e-06
256 4.1963e-06 4.0520e-06
Table 3: The time integral error between the 4-order integral scheme and the Chebyshev integral method

4.1.2 The efficiency of SCFT calculations

To further demonstrate the performance of our proposed approach, we apply the numerical schemes to SCFT calculations. To compare results, we need a metric for accuracy that can be readily compared across different calculations. We use the value of the single chain partition function, , as the metric of accuracy of the solver. Since it is an integration of the end result of the propagator solve, it is a measure of the entire solution process. As a basis for comparison, we use a square with the edge length of as the computational domain. The volume fraction of is , and the interaction parameter . The computation is carried out using a quadrilateral mesh (see Fig. 4 (a)). Correspondingly, the convergent morphology is a cylindrical structure, as shown in Fig. 4 (b).

Figure 4: Cylindrical phase calculated by VEM with uniform grid when , . Red colors correspond to large A-segment fractions.

First, we look at the contour discretization schemes. The goal is to have the fewest number of contour points necessary for a desired accuracy. The quadratic VEM with nodes is used to guarantee enough spatial discretization accuracy. in Fig. 5(a) is numerically obtained by the SDC scheme with contour points. Fig. 5(a) shows the convergence of for the CN and SDC schemes as discussed above. The SDC method converges faster than the CN scheme to a prescribed precision.

Second we observe the numerical behavior of linear and quadratic VEMs in the SCFT simulation. From the above numerical tests (see Fig. 5(a)), one can see that using SDC scheme with discretization points can guarantee enough accuracy in time direction. So in the following computations, we use a high-precision numerical as the exact value, which is obtained by the quadratic VEM with nodes and SDC scheme with points. Fig. 5 (b) shows the values with different spatial discretization points of linear and quadratic VEMs. It is easy to see that the quadratic VEM is more accurate than the linear VEM as theory predicts. Therefore, in the following calculations, we always adopt the quadratic VEM and the SDC scheme.

Figure 5: Convergence of numerically computing single chain partition function by different schemes. is the relative error. is the numerical exact solution . (See text for the details about ). (a) shows results obtained by the CN and SDC scheme as the contour points increase. Quadratic VEM with is employed to discretize the spatial variable. (b) presents results computed by the linear and quadratic VEMs with a increase of spatial discretization points. SDC scheme with nodes is applied to discretize the time variable.

4.2 General domains with general polygonal meshes

One advantage of the VEM is the capacity of approximating arbitrary geometry domain with general polygonal meshes. This allows us to directly calculate ordered structures on physical domains. Fig. 6 presents these results on five different two dimensional domains discreted by quadranglar and polygonal elements, respectively. The same convergent structure and almost the same Hamiltonian value can be obtained for these two kinds of meshes, as shown in Fig.6 and in Tab. 4.

Figure 6: The self-assembled patterns in general domains through SCFT simulation including (1). Flower shaped plane; (2). Curved-L shaped plane; (3). Ring domain; (4). Rabbit-shaped plane; and (5). Dumbbell plane. Red colors correspond to large A-segment fractions. The first and second columns present the schematic mesh of quadrangular and polygonal meshes, respectively. The simulating diblock copolymer systems contain (1c) , (1d) , (2c) , (2d) , (3c) , (3d) , (4c) , (4d) , (5c) , and (5d) . The number of used nodes and converged Hamiltonian values can be found in Tab. 4.
Patterns Grid Nodes Hamiltonian
(c) (d)
(1) (a) 13041 -2.3742 -1.7388
(b) 22560 -2.3754 -1.7398
(2) (a) 10720 -2.3720 -1.7440
(b) 20273 -2.3765 -1.7382
(3) (a) 7014 -3.1410 -0.1874
(b) 6510 -3.1409 -0.1873
(4) (a) 30182 -3.1440 -0.1900
(b) 34587 -3.1448 -0.1901
(5) (a) 7601 -2.3670 -1.6797
(b) 13824 -2.3718 -1.6883
Table 4: The number of nodes of different meshes used in SCFT calculations for five different domains as shown in Fig. 6 and corresponding converged Hamiltonian values.

4.3 Adaptive VEM

In this subsection, we will demonstrate the efficiency of adaptive VEM from two parts: 1) spending less computational cost to obtain prescribed accuracy; 2) application to strong segregation systems. As discussed in Sec. 4.1 the quadratic VEM is more accurate than the linear one, therefore, only the quadratic VEM is used in the adaptive process. Meanwhile the SDC scheme with contour discretization points is chosen in time direction.

First, we take and as an example to demonstrate the efficiency of the adaptive method. The computational domain is a square with edge length . The uniform mesh of square cells with nodes is used to model the system at the start stage, then the adaptive method is launched when the iteration reaches the maximum steps or the reference value of estimator (,

is the standard deviation of

, estimator see Eqn. (56)). The adaptive process will be terminated when the successive Hamiltonian difference is smaller than . Fig. 7 (a) gives the final adaptive meshes which includes nodes. Fig. 7 (b) shows the convergent tendency of Hamiltonian of adaptive process. The finally converged morphology has been shown in Fig. 4 (b). It can be seen that the Hamiltonian value efficiently converges by the cascadic adaptive grid method and grid elements adaptively increase mainly concentrated on the shape interface.

Figure 7: (a) The converged adaptive meshes. (b) The numerical behavior of Hamiltonian . The numbers between two dotted lines represent the number of spatial nodes in the adaptive process.

We also compared the simulation results on the adaptive mesh and the uniform mesh. Fig. 8 shows the numerical behaviors of the single chain partition function and Hamiltonian as the number of mesh nodes increase. Tab. 5 gives the corresponding converged values of and on the adaptive mesh and the uniform mesh with nodes, respectively. From these results, one can find that the uniform mesh results indeed gradually converge to the adaptive results. However there exists a small gap between the results of adaptive mesh and that of uniform mesh. The reason is that, compared with the uniform mesh method, the adaptive method put more mesh nodes in the areas where the solution changes sharply. The minimum element size of the adaptive mesh in the above calculation is . If the same element size is used, the uniform mesh requires about nodes which is about ten times than the adaptive method.

Figure 8: The convergence results on the adaptive mesh and the uniform mesh when . The differences of (a) single partition function , and (b) Hamiltonian value . and are obtained on the uniform mesh, while and are adaptive results.
Mesh Nodes Q H
Adaptive 6684 4.2295e+02 -2.369403
Uniform 16641 4.2373e+02 -2.369448
Table 5: The values of and .

Next, we apply the adaptive method to simulate the strong segregation systems, i.e., large interaction parameter , also in the square domain with edge length . For strong segregation case, the interface thickness becomes narrower. Therefore the adaptive mesh is more suitable than the uniform mesh to catch these narrower interface. When simulating the strong segregation system, the initial values are obtained by the converged results of relatively weak segregation system. Tab. 6 presents the numerical results of from to and . From these results, one can find the advantages of the adaptive method as increases, including a mild increase of mesh nodes, and a less iteration steps. Finally, we apply the adaptive method to strong segregation systems on more complicated domains, including two kinds of structures, cylindrical structures when , and lamellar phases when , . Fig. 9 presents the adaptive meshes and converged morphologies. These results demonstrate that the adaptive VEM method indeed can improve precision and reduce the computational cost.

Step Nodes Q H
25 1146 6684 4.2295e+02 -2.369403
30 78 9037 3.7743e+03 -3.149607
35 89 13443 3.1519e+04 -4.020791
40 74 17649 2.5880e+05 -4.946249
45 75 19741 2.0408e+06 -5.907039
50 75 20480 1.5355e+07 -6.892386
55 73 20641 1.0513e+08 -7.895548
60 61 20690 6.5573e+08 -8.911902
Table 6: Numerical results by the adaptive VEM for strong segregation systems.
Figure 9: The self-assembled patterns of strong segregation systems obtained by the adaptive VEM. Red colors correspond to large A-segment fractions, including (1b) , (1d) , (2b) , and (2d) . The first and third columns present the corresponding adaptive meshes.

5 Conclusion

In this paper, we propose a high-accuracy numerical method to solve the polymer SCFT model on arbitrary domains. The approach combines high order VEM and SDC scheme together for solving the propagator equations, which is the most-consuming part of the SCFT simulation. The VEM can use on very general polygonal mesh and is more flexible to approximate complex computational domains. We also develop an adaptive equipped with a new marking strategy to efficiently make use of the information of existing numerical results and to significantly save the SCFT iterations. The resulting method achieves high accuracy with fewer number of spatial and contour nodes, and is suitable for solving strong segregation systems. Numerical results demonstrate the efficiency of our proposed approach.


This work is supported by National Natural Science Foundation of China (11771368, 11871413), Hunan Science Foundation of China (2018JJ2376), Youth Project (18B057) and Key Project (19A500) of Education Department of Hunan Province of China.


  • (1) C. R. Stewart-Sloan and E. L. Thomas. Interplay of symmetries of block polymers and confining geometries. Eur. Polym. J. 47 (04) (2011) 630-646.
  • (2) R. A. Segalman. Patterning with block copolymer thin films. Mater. Sci. Eng., R 46 (06) (2005) 191-226.
  • (3) H. Wang, B. Shentu, R. Faller. Molecular dynamics of different polymer blends containing poly(2,6-dimethyl-1,4-phenylene ether). Phys. Chem. Chem. Phys. 17 (06) (2015) 4714-4723.
  • (4) V. Sethuraman, B.H. Nguyen, V. Ganesan. Coarse-graining in simulations of multicomponent polymer systems. J. Chem. Phys. 141 (2014) 244904.
  • (5) G. H. Fredrickson. The equilibrium theory of inhomogeneous polymers. Oxford University Press: New York, (2006).
  • (6) W. Xu, K. Jiang, P. Zhang, A. C. Shi. A strategy to explore stable and metastable ordered phases of block copolymers. J. Phys. Chem. B 117 (17) (2013) 5296-5305.
  • (7) K. Jiang, Y. Huang, P. Zhang. Spectral method for exploring patterns of diblock copolymers. J. Comput. Phys. 229 (20) (2010) 7796-7805.
  • (8) K. Jiang, C. Wang, Y. Huang, P. Zhang. Discovery of new metastable patterns in diblock copolymers. Commun. Comput. Phys. 14 (02) (2013) 443-460.
  • (9) G. Ouaknin, N. Laachi, K. Delaney, G. H. Fredrickson. Self-consistent field theory simulations of polymers on arbitrary domains. J. Comput. Phys. 327 (2016) 168-185.
  • (10) H. D. Ceniceros. Efficient order-adaptive methods for polymer self-consistent field theory. J. Comput. Phys. 386 (01) (2019) 9-21.
  • (11) H. D. Ceniceros, G. H. Fredrickson. Numerical solution of polymer self-consistent field theory. Multiscale Model. Simul. 2 (03) (2004) 452-474.
  • (12) R. B. Thompson, K. Ø. Rasmussen, T. Lookman. Improved convergence in block copolymer self-consistent field theory by anderson mixing. J. Chem. Phys. 120 (31) (2004) 31-34.
  • (13) K. Jiang, W. Xu, P. Zhang. Analytic structure of the SCFT energy functional of multicomponent block copolymers. Commun. Comput. Phys. 17 (05) (2015) 1360-1387.
  • (14) M. W. Matsen, M. Schick. Stable and unstable phases of a diblock copolymer melt. Phys. Rev. Lett. 72 (16) (1994) 2660-2663.
  • (15) K. Ø. Rasmussen, G. Kalosakas. Improved numerical algorithm for exploring block copolymer mesophases. J. Phys.: Condens. Matter 40 (16) (2002) 1777-1783.
  • (16) E. W. Cochran, C.J. Garcia-Cervera, G. H. Fredrickson. Stability of the gyroid phase in diblock copolymers at strong segregation. Macromolecules 39 (07) (2006) 2449-2451.
  • (17) D. M. Ackerman, K. Delaney, G. H. Fredrickson, B. Ganapathysubramaniana. A finite element approach to self-consistent field theory calculations of multiblock polymers. J. Comput. Phys. 331 (2017) 280-296.
  • (18) H. Wei, M. Xu, W. Si, and K. Jiang. A finite element method of the self-consistent field theory on general curved surfaces. J. Comput. Phys. 387 (15) (2019) 230-244.
  • (19) P. Binev, W. Dahmen, and R. DeVore. Adaptive finite element methods with convergence rates. Numer. Math. (Heidelb), 97 (02) (2004) 219-268.
  • (20) C. Carstensen, J. Hu. Hanging nodes in the unifying theory of a posteriori finite element error control. J. Comput. Math., (2009) 215-236.
  • (21) B. Ahmed, A. Alsaedi, F. Brezzi, L.D. Marini, and A. Russo. Equivalent projectors for virtual element methods. Comput. Math. Appl., 66 (03) (2013) 376-391.
  • (22) L. Beirao da Veiga, F. Brezzi, L.D. Marini, and A. Russo. The hitchhiker’s guide to the virtual element method. Math. Models Methods Appl. Sci. 24 (08) (2014) 1541-1573.
  • (23) L. Beirao da Veiga, F. Brezzi, A. Cangiani, L. D. Marini, G. Manzini and A. Russo. Basic principles of virtual elements methods. Math. Models Methods Appl. Sci. 23 (2013) 199-214.
  • (24) P. F. Antonietti, L. Beirão da Veiga, D. Mora, and M. Verani. A stream function formulation of the Stokes problem for the virtual element method. SIAM J. Numer. Anal. 52 (2014) 386–404.
  • (25) J. Zhao, S. Chen, and B. Zhang. The nonconforming virtual element method for plate bending problems. Math. Models Methods Appl. Sci. 26 (2016) 1671–1687.
  • (26) L. Chen, H. Wei, and M. Wen. An interface-fitted mesh generator and virtual element methods for elliptic interface problems. J. Comput. Phys. 334 (2017) 327–348.
  • (27)

    A. Dutt, L. Greengard, and V. Rokhlin. Spectral deferred correction methods for ordinary differential equations. BIT Numer. Math., 40 (02) (2000) 241-266.

  • (28) C. W. Clenshaw and A. R. Curtis. A method for numerical integration on an automatic computer. Numer. Math. (Heidelb), 2 (01) (1960) 197-205.
  • (29) U. M. Ascher, S. J. Ruuth, and R. J. Spiteri. Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Appl. Numer. Math., 25 (1997) 151-167.
  • (30) L. N. Trefethen. Is Gauss quadrature better than Clenshaw-Curtis? SIAM Review, 50 (01) (2008) 67-87.
  • (31) Y. Huang, K. Jiang, N. Yi. Some weighted averaging methods for gradient recovery. Adv. Appl. Math. Mech. 4 (02) (2012) 131-155.
  • (32) H. Jarausch. On an adaptive grid refining technique for finite element approximations. SIAM J. Sci. and Stat. Comput., 7(04) (1986) 1105-1120.
  • (33) W. Dörfler. A convergent adaptive algorithm for poisson’s equation. SIAM J. Numer. Anal., 33(03) (1996) 1106-1124.
  • (34) H. Wei and K. Jiang. A new efficient marking strategy for adaptive methods, in preparation.
  • (35) S. J. Osher, R. P. Ronald. Level set methods and dynamic implicit surfaces. New York: Springer, (2005).
  • (36) H. Wei. FEALPy: finite element analysis library in Python. https://github.com/weihuayi/