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, confinementinduced entropy loss, and surface interactions, can lead to novel morphologies which are not obtained in the bulk. There are many industrial applications for the blockcopolymer ordered structures at the nanoscale, such as the construction of highcapacity 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 coarsegraining 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 selfconsistent field theory (SCFT), which is one of the most successful modern tools for studying the phase behaviors of inhomogeneous polymer systems, such as selfassembly, 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 particlebased model into a fieldtheoretic framework, finally obtains a meanfield equations system within saddlepoint approximation fredrickson2006equilibrium .
From the viewpoint of mathematics, the SCFT model is a complicated variational problem with many challenges, such as saddlepoint, nonlinearity, multisolutions, and multiparameters. 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 timedependent partial differential equations (PDEs)
ouaknin2016scft ; ceniceros2019efficient , evaluating (monomer) density operators ceniceros2019efficient , and finding saddlepoints 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 highorder 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 Selfconsistent 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 fieldbased Hamiltonian within meanfield approximation for the incompressible diblock copolymer melt is fredrickson2006equilibrium ; cochran2006stability
(1) 
where is the FloryHuggins 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
(2) 
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(3a)  
(3b) 
with the initial condition and being the radius of gyration. The above PDE is welldefined through possessing an appropriate boundary condition. In this work, we consider the homogeneous Neumann boundary condition
(4) 
The reverse propagator , which represents the probability weight from to , satisfies Eqn. (3) only with the righthand 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
(5)  
(6) 
The first variations of the Hamiltonian with respect to fields and lead to the meanfield equations
(7)  
(8) 
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 selfconsistent iterative procedure as shown in the following flowchart.
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, highaccuracy 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
(9)  
An accelerated semiimplicit scheme has been developed to find the equilibrium states ceniceros2004numerical ; jiang2015analytic
. However, the existing semiimplicit 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 timeconsuming 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 nonconvex) 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.
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 nonoverlapping 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
(10) 
and . We will also use instead of , where is a onedimensional index of the nature correspondence of , for example,
(11) 
The local virtual element space can be defined as B2013vem ; lfal2013vem
(12) 
where denotes the common Laplace operator. is a set of polynomials of degree up to on . The dimension of is
(13) 
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 GaussLobatto quadrature points on e;
Then the global virtual element space can be defined based on the local space
(14) 
The dimension of is
(15) 
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
(16) 
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 ,
(17) 
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
(18) 
Let , using the expression (16), . The discretized variational formulation (18) has the matrix form
(19) 
where
and
(20) 
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 ,
(21) 
which projects the local virtual element space onto the polynomial space with degree up to . For each , we have the orthogonality condition
(22) 
The above condition defines only up to a constant. It can be fixed by prescribing a projection operator onto constants requiring
(23) 
can be chosen as
(24)  
where is the number of vertices of .
Next we compute local stiffness matrix on the polygon ,
(25) 
With the projector , can be splitted into
(26) 
then Eqn. (25) becomes
(27) 
Replacing the second term as
we can obtain the approximate local stiffness matrix
(28) 
3.1.4 Construction of the mass matrix
The mass matrix in the VEM can be obtained from the local projection, . For each ,
(29) 
can not be calculated directly. Next, we show how to compute the local mass matrix L2013vem
(30) 
Similar to the construction method of stiffness matrix, we can define the basis function through projection operator
(31) 
Then
(32) 
Replacing the second term in the above equation as
the local mass matrix can be approximated as
(33) 
3.1.5 Cross mass matrix
The local cross mass matrix on can be defined as
(34) 
Applying the projection , as defined in the above Sec. 3.1.4, into the cross term, the local mass matrix can be calculated as
(35) 
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.
(36) 
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 Chebyshevnode interpolatory quadrature, to integrate the contour derivative, which can achieve a highaccuracy 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 CrankNicholson (CN) scheme to solve the semidiscretize propagator equation (19), and obtain the initial numerical solution .
(37) 
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 semidiscretize propagator equation (19), such as secondorder operatorsplitting method rasmussen2002improved , implicitexplicit RungeKutta scheme ascher1997rungekutta .
Then we use the deferred correction scheme to achieve a highaccuracy numerical solution. We can give the exact semidiscretize solution of propagator by integrating (19) along the contour variable
(38) 
The error between the numerical solution and exact semidiscretize solution is defined as
(39) 
Multiplying both sides by , we have
(40)  
(41)  
(42)  
(43) 
the residual
(44) 
can be compute by the spectral integral method with Chebyshevnodes as presented in the Sec. 3.3. By the definition of residual , we have the error integration equation
(45) 
Taking the first derivative of the above equation with respect to leads to
(46) 
which can be also solved by the CN scheme (37). Then the corrected numerical solution is
(47) 
Repeating the above process, one can have , is the predetermined number of deferred corrections. The convergent order of deferred correction solution along the contour parameter is
(48) 
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
,(49) 
is the true solution of propagator, and .
3.3 Spectral integral method along the contour variable
In this section, we discuss the Chebyshevnode 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
(50) 
the integrand is a smooth function. After changing variables, the general integral becomes
(51) 
where
. The interpolate polynomial of
at Chebyshev nodes , .(52) 
where
(53) 
. In practice, coefficients are calculated by the fast discrete cosine transform.
(54) 
Due to , we have
(55) 
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 ,
(56) 
is the harmonic average operator huang2012
(57) 
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.
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
(58) 
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.
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
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.7737e02  –  1.2062e03  – 
1089  1.3267e02  1.84  1.5084e04  2.99 
4225  3.4013e03  1.96  1.8863e05  2.99 
16641  8.5563e04  1.99  2.3582e06  3.00 
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 .
CN  SDC  
order  order  
4  6.0605e03  –  5.7514e04  – 
8  1.5074e03  2.00  1.0163e05  5.82 
16  3.7637e04  2.00  6.4626e07  3.97 
32  9.4065e05  2.00  4.2283e08  3.94 
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 semidiscrete system. Correspondingly the exact solution of (59) can be discretized into . Then we solve the semidiscretize 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
(60) 
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.2343e03  2.4219e04 
8  1.0728e03  4.4032e06 
16  2.6970e04  4.0527e06 
32  6.7675e05  4.0519e06 
64  1.7401e05  4.0520e06 
128  5.8778e06  4.0520e06 
256  4.1963e06  4.0520e06 
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).
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 highprecision 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.
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.
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 
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.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.
Mesh  Nodes  Q  H 
Adaptive  6684  4.2295e+02  2.369403 
Uniform  16641  4.2373e+02  2.369448 
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 
5 Conclusion
In this paper, we propose a highaccuracy 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 mostconsuming 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.
Acknowledgements
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.
References
 (1) C. R. StewartSloan and E. L. Thomas. Interplay of symmetries of block polymers and confining geometries. Eur. Polym. J. 47 (04) (2011) 630646.
 (2) R. A. Segalman. Patterning with block copolymer thin films. Mater. Sci. Eng., R 46 (06) (2005) 191226.
 (3) H. Wang, B. Shentu, R. Faller. Molecular dynamics of different polymer blends containing poly(2,6dimethyl1,4phenylene ether). Phys. Chem. Chem. Phys. 17 (06) (2015) 47144723.
 (4) V. Sethuraman, B.H. Nguyen, V. Ganesan. Coarsegraining 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) 52965305.
 (7) K. Jiang, Y. Huang, P. Zhang. Spectral method for exploring patterns of diblock copolymers. J. Comput. Phys. 229 (20) (2010) 77967805.
 (8) K. Jiang, C. Wang, Y. Huang, P. Zhang. Discovery of new metastable patterns in diblock copolymers. Commun. Comput. Phys. 14 (02) (2013) 443460.
 (9) G. Ouaknin, N. Laachi, K. Delaney, G. H. Fredrickson. Selfconsistent field theory simulations of polymers on arbitrary domains. J. Comput. Phys. 327 (2016) 168185.
 (10) H. D. Ceniceros. Efficient orderadaptive methods for polymer selfconsistent field theory. J. Comput. Phys. 386 (01) (2019) 921.
 (11) H. D. Ceniceros, G. H. Fredrickson. Numerical solution of polymer selfconsistent field theory. Multiscale Model. Simul. 2 (03) (2004) 452474.
 (12) R. B. Thompson, K. Ø. Rasmussen, T. Lookman. Improved convergence in block copolymer selfconsistent field theory by anderson mixing. J. Chem. Phys. 120 (31) (2004) 3134.
 (13) K. Jiang, W. Xu, P. Zhang. Analytic structure of the SCFT energy functional of multicomponent block copolymers. Commun. Comput. Phys. 17 (05) (2015) 13601387.
 (14) M. W. Matsen, M. Schick. Stable and unstable phases of a diblock copolymer melt. Phys. Rev. Lett. 72 (16) (1994) 26602663.
 (15) K. Ø. Rasmussen, G. Kalosakas. Improved numerical algorithm for exploring block copolymer mesophases. J. Phys.: Condens. Matter 40 (16) (2002) 17771783.
 (16) E. W. Cochran, C.J. GarciaCervera, G. H. Fredrickson. Stability of the gyroid phase in diblock copolymers at strong segregation. Macromolecules 39 (07) (2006) 24492451.
 (17) D. M. Ackerman, K. Delaney, G. H. Fredrickson, B. Ganapathysubramaniana. A finite element approach to selfconsistent field theory calculations of multiblock polymers. J. Comput. Phys. 331 (2017) 280296.
 (18) H. Wei, M. Xu, W. Si, and K. Jiang. A finite element method of the selfconsistent field theory on general curved surfaces. J. Comput. Phys. 387 (15) (2019) 230244.
 (19) P. Binev, W. Dahmen, and R. DeVore. Adaptive finite element methods with convergence rates. Numer. Math. (Heidelb), 97 (02) (2004) 219268.
 (20) C. Carstensen, J. Hu. Hanging nodes in the unifying theory of a posteriori finite element error control. J. Comput. Math., (2009) 215236.
 (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) 376391.
 (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) 15411573.
 (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) 199214.
 (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 interfacefitted 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) 241266.
 (28) C. W. Clenshaw and A. R. Curtis. A method for numerical integration on an automatic computer. Numer. Math. (Heidelb), 2 (01) (1960) 197205.
 (29) U. M. Ascher, S. J. Ruuth, and R. J. Spiteri. Implicitexplicit RungeKutta methods for timedependent partial differential equations. Appl. Numer. Math., 25 (1997) 151167.
 (30) L. N. Trefethen. Is Gauss quadrature better than ClenshawCurtis? SIAM Review, 50 (01) (2008) 6787.
 (31) Y. Huang, K. Jiang, N. Yi. Some weighted averaging methods for gradient recovery. Adv. Appl. Math. Mech. 4 (02) (2012) 131155.
 (32) H. Jarausch. On an adaptive grid refining technique for finite element approximations. SIAM J. Sci. and Stat. Comput., 7(04) (1986) 11051120.
 (33) W. Dörfler. A convergent adaptive algorithm for poisson’s equation. SIAM J. Numer. Anal., 33(03) (1996) 11061124.
 (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/