Goal-Oriented Adaptive Mesh Refinement for the Quasicontinuum Approximation of a Frenkel-Kontorova Model
aa r X i v : . [ m a t h . NA ] N ov GOAL-ORIENTED ADAPTIVE MESH REFINEMENT FOR THEQUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVAMODEL
MARCEL ARNDT AND MITCHELL LUSKIN
Abstract.
The quasicontinuum approximation [24] is a method to reduce the atomistic degrees offreedom of a crystalline solid by piecewise linear interpolation from representative atoms that arenodes for a finite element triangulation. In regions of the crystal with a highly nonuniform deforma-tion such as around defects, every atom must be a representative atom to obtain sufficient accuracy,but the mesh can be coarsened away from such regions to remove atomistic degrees of freedom whileretaining sufficient accuracy. We present an error estimator and a related adaptive mesh refinementalgorithm for the quasicontinuum approximation of a generalized Frenkel-Kontorova model thatenables a quantity of interest to be efficiently computed to a predetermined accuracy. Introduction
The solution of the equations for mechanical equilibria of a crystalline solid modeled by a classicalatomistic potential requires the computation of the interaction of each atom with all of the otheratoms in its sphere of influence. Due to the high computational complexity, it is generally notpossible to obtain numerical solutions for systems that are large enough to simulate long-rangeelastic effects, even for short-ranged potentials. However, the local environment of nearby atomsis almost identical up to translation, except in the neighborhood of defects such as cracks anddislocations. The quasicontinuum method utilizes this slow variation of the strain away fromdefects to approximate the full systems of equations of mechanical equilibrium by equations ofequilibrium at a reduced set of representative atoms [7, 8, 16, 24].More precisely, the positions of the full set of atoms are obtained by piecewise linear interpolationfrom representative atoms that are nodes for a finite element triangulation. A quasicontinuumenergy is defined as a function of the positions of the representative atoms. The quasicontinuummethod is then used to obtain a solution to a desired accuracy with a significant reduction in thecomputational degrees of freedom by coarsening the finite element mesh away from the singularities.Near defects, sufficient accuracy can only be obtained if all atoms are representative atoms. Incontrast to continuum models, the mesh cannot be refined past the atomistic scale.Even higher efficiency is achieved by approximating the total energy of the atoms in a coarsenedtriangle by the product of the strain energy density and the area of the triangle (or its higherdimensional analogue). The strain energy density is obtained from the energy per atom in a latticewhich is a uniform strain of the infinite reference lattice. The uniform strain in turn is determinedfrom the displacement of the nodes of the respective triangle. The quasicontinuum approximationwe obtain this way allows the coupling between a region that is computed as in fully atomistic
Mathematics Subject Classification.
Key words and phrases.
Error estimation, a posteriori , adaptive, refinement, goal-oriented, atomistic-continuummodeling, quasicontinuum, Frenkel-Kontorova model, dislocation, defect.This work was supported in part by DMS-0304326 and by the Minnesota Supercomputing Institute. This work isalso based on work supported by the Department of Energy under Award Number DE-FG02-05ER25706. simulations and a region that is computed using the methods of piecewise linear finite elementcontinuum mechanics.It is well-known that an efficient adaptive algorithm is highly dependent on the quantity ofinterest or goal of the computation, see [1, 5]. In this paper, we adopt this goal-oriented approachto obtain the approximation of a quantity of interest to within a desired tolerance.The reliable and efficient utilization of the quasicontinuum method requires both a strategy todetermine the decomposition of the fully atomistic system into atomistic and continuum regionsand a strategy for refinement within the continuum region. In [2, 3], we have developed a goal-oriented error estimator and a corresponding adaptive algorithm to decide between the atomisticmodel and the continuum model. In this paper, we develop an error estimator and a correspondingadaptive algorithm for mesh refinement in the continuum region.We have chosen initially to investigate as a model problem a generalization of the classicalFrenkel-Kontorova model. The potential energy includes next-nearest-neighbor interactions in ad-dition to nearest-neighbor interactions so that the continuum energy of a representative atom isdifferent from the atomistic energy of a representative atom for a nonuniform strain.We note that in contrast to error estimators and adaptive algorithms for mesh refinement inclassical continuum mechanics, our continuum region is coupled to the atomistic region. This isa considerably more complex setting than in purely continuum models with classical boundaryconditions. Also, our mesh refinement algorithm restricts the representative atoms which serve asthe mesh points to the sites of the atoms in the reference lattice, although this could be relaxedaway from the atomistic-continuum interfacial region.Algorithms for adaptive mesh refinement for the quasicontinuum method have been proposed andinvestigated in numerical experiments for several mechanics problems in [11, 17, 23]. An a posteriori error indicator for a global norm was analyzed and tested for a variant of our one-dimensionalquasicontinuum method in [20]. The goal-oriented approach to adaptive mesh refinement for thequasicontinuum method was first investigated in [18, 19]. Mathematical analyses of several variantsof the quasicontinuum method have been given in [6, 9, 10, 12, 13, 21]. We refer to [4, 7, 14, 22, 25]for alternative atomistic-continuum coupling methods.2.
Quasicontinuum Approximation
In this section, we introduce the atomistic model and its quasicontinuum approximation. Werefer to [2] for a more detailed description.We consider a one-dimensional system of 2 M atoms whose positions are denoted by y =( y − M +1 , . . . , y M ) ∈ R M . The potential energy of the atomistic system is described by a func-tion E a : R M → R . (2.1)We split the energy into atom-wise contributions by means of E a ( y ) = M X i = − M +1 E ai ( y ) . (2.2)In this paper, we consider a Frenkel-Kontorova type model [15] which serves as a one-dimensionaldescription of a dislocation. We expect that the a posteriori error estimators that we introduce andthe corresponding adaptive refinement algorithms will be applicable to more general quasicontinuummodels. For the Frenkel-Kontorova model, the atom-wise contribution, E ai , consists of two parts, E ai = E a,ei + E a,mi . (2.3) DAPTIVE REFINEMENT FOR THE QUASICONTINUUM APPROXIMATION 3 −2a −3a −4a −5a −6a −a substratemodeledlayermisfitenergy Figure 1:
Frenkel-Kontorova model. The wells depict the misfit energy, E m,ai . The elastic part, E a,ei , describes nearest-neighbor (NN) interactions and next-nearest-neighbor(NNN) interactions, whereas E a,mi models the misfit energy of a slip plane sitting on an unde-formed substrate, see Figure 1. The two parts are defined as E a,ei ( y ) = k ( y i − y i − − a ) + k ( y i +1 − y i − a ) + k ( y i − y i − − a ) + k ( y i +2 − y i − a ) , E a,mi ( y ) = ( k ( y i − ( i − a ) , i = − M + 1 , . . . , , k ( y i − ia ) , i = 1 , . . . , M. (2.4)Here a ∈ R denotes the equilibrium distance, and the moduli k , k , and k describe the strength ofthe misfit energy and the elastic interactions, respectively. To ensure coercivity, we require k > k + 2 k > | k | .Next, the continuum energy E ci = E c,ei + E c,mi (2.5)for each atom i is derived from the atomistic energy, which leads to E c,ei ( y ) = k ( y i − y i − − a ) + k ( y i +1 − y i − a ) , E c,mi ( y ) = ( k ( y i − ( i − a ) , i = − M + 1 , . . . , , k ( y i − ia ) , i = 1 , . . . , M, (2.6)where k = k + 4 k , see [2] for the details. We note that E c,ei ( y ) = E a,ei ( y ) if y i +2 − y i +1 = y i +1 − y i = y i − y i − = y i − − y i − , but E c,ei ( y ) = E a,ei ( y ) more generally.For each atom i , we decide whether this atom is modeled atomistically or as continuum. An aposteriori error estimator for this task has been derived in [2]. Let δ ai = ( i is modeled atomistically,0 if atom i is modeled as continuum, and δ ci = 1 − δ ai . (2.7)We define the atomistic-continuum energy to be E ac ( y ) := M X i = − M +1 δ ai E ai ( y ) + M X i = − M +1 δ ci E ci ( y ) . (2.8) DAPTIVE REFINEMENT FOR THE QUASICONTINUUM APPROXIMATION 4
The next step is to coarsen out unnecessary atoms within the continuum region. This way, werestrict the system to the remaining atoms, called the representative atoms, or briefly repatoms.An a posteriori error estimator to determine the optimal coarsening will be developed in this paper.Let ℓ j for j = − N +1 , . . . , N be the index of the j -th repatom, where 2 N is the number of repatoms.We require that ℓ − N +1 < ℓ − N +2 < · · · < ℓ j < ℓ j +1 < · · · < ℓ N − < ℓ N (2.9)and ℓ − N +1 = − M + 1 , ℓ − N +2 = − M + 2 , ℓ N − = M − , ℓ N = M. (2.10)Then ν j := ℓ j +1 − ℓ j (2.11)gives the number of atomistic intervals between the repatoms j and j + 1. We denote the vector ofrepatoms by y qc .The quasicontinuum energy is obtained by implicitly reconstructing the missing atoms from therepatoms by piecewise linear interpolation, and then computing the atomistic-continuum energyfrom the reconstructed vector. The piecewise linear interpolation can be written as the matrixmultiplication I aq y qc , where I aqℓ j + k,j := ν j − kν j , I aqℓ j + k,j +1 := kν j (2.12)for j = − N + 1 , . . . , N and k = 0 , . . . , ν j , and I aqi,j = 0 otherwise. Hence, the quasicontinuum energyis given by E qc ( y qc ) := E ac ( I aq y qc ) . (2.13)Summation formulas for an efficient computation of E qc ( y qc ) without having to explicitly recon-struct the non-repatoms have been derived in [2].We describe boundary conditions at the two leftmost atoms and the two rightmost atoms. Wetake into account two instead of one atom at each end of the chain because of the NNN interactions.For the fully atomistic system and the atomistic-continuum system, the boundary conditions readas y a − M +1 = y bcl , y a − M +2 = y bcl , y aM − = y bcr , y aM = y bcr ,y ac − M +1 = y bcl , y ac − M +2 = y bcl , y acM − = y bcr , y acM = y bcr , (2.14)for given values y bcl , y bcl , y bcr , and y bcr , and similarly the boundary conditions for the quasicontinuumsystem are given by y qc − N +1 = y bcl , y qc − N +2 = y bcl , y qcN − = y bcr , y qcN = y bcr . (2.15)Next, we introduce the solution spaces V a := R M , V a := R M − , V q := R N , V q := R N − (2.16)and the boundary operators J a : V a → V a and J q : V q → V q (2.17)by J ai,j := δ ij for i = − M + 1 , . . . , M, j = − M + 3 , . . . , M − ,J qi,j := δ ij for i = − N + 1 , . . . , N, j = − N + 3 , . . . , N − . (2.18) DAPTIVE REFINEMENT FOR THE QUASICONTINUUM APPROXIMATION 5 V q V a V q V a / / I aq O O J q O O J a Figure 2:
Solution spaces and their operators for the fully-atomistic level (a) and the quasicontinuum level(q). V a and V a are the suitable spaces for the atomistic system and the atomistic-continuum system,both with and without boundary values, whereas V q and V q are the suitable spaces for the qua-sicontinuum system. J a and J q extend vectors from V a and V q , respectively, by zero boundaryvalues. Note that the interpolation operator I aq maps V q to V a . The spaces and their operatorsare illustrated in Figure 2.To implement the boundary conditions, we define the vectors y bcq := (cid:2) y bcl y bcl · · · y bcr y bcr (cid:3) T ∈ V q , y bca := I aq y bcq ∈ V a . (2.19)We seek minimizers of the three different potential energies subject to the boundary conditionsdescribed above, that is, vectors ¯y a ∈ J a V a + y bca ⊂ V a , ¯y ac ∈ J a V a + y bca ⊂ V a , and ¯y qc ∈ J q V q + y bcq ⊂ V q (2.20)which minimize the potential energies E a , E ac , and E qc , respectively. If we decompose ¯y a = J a y a + y bca , ¯y ac = J a y ac + y bca , and ¯y qc = J q y qc + y bcq , (2.21)then the fully atomistic solution y a ∈ V a , the atomistic-continuum solution y ac ∈ V a , and thequasicontinuum solution y qc ∈ V q are characterized by y a = arg min y ∈ V a E a ( J a y + y bca ) , y ac = arg min y ∈ V a E ac ( J a y + y bca ) , y qc = arg min y ∈ V q E qc ( J q y + y bcq ) . (2.22)Next, we write the energies in matrix notation as E a ( y ) = ( y − a a ) T D T E a D ( y − a a ) + ( y − b a ) T K ( y − b a ) , E ac ( y ) = ( y − a a ) T D T E ac D ( y − a a ) + ( y − b a ) T K ( y − b a ) , E qc ( y ) = ( I aq y − a a ) T D T E ac D ( I aq y − a a ) + ( I aq y − b a ) T K ( I aq y − b a ) . (2.23)We refer to [2] for the precise and lengthy definitions of the respective matrices. In matrix notation,the vectors y a , y ac , and y qc are the solutions of the linear systems M a y a = f a ,M ac y ac = f ac ,M qc y qc = f qc , (2.24) DAPTIVE REFINEMENT FOR THE QUASICONTINUUM APPROXIMATION 6 where the matrices M a , M ac , and M qc are given by M a := J a T ( D T E a D + K ) J a ,M ac := J a T ( D T E ac D + K ) J a ,M qc := J q T I aq T ( D T E ac D + K ) I aq J q , (2.25)and where the right-hand sides f a , f ac , and f qc are defined as f a := − J a T D T E a D ( y bca − a a ) − J a T K ( y bca − b a ) , f ac := − J a T D T E ac D ( y bca − a a ) − J a T K ( y bca − b a ) , f qc := − J q T I aq T D T E ac D ( I aq y bcq − a a ) − J q T I aq T K ( I aq y bcq − b a ) . (2.26)3. Goal-Oriented Error Estimation
We estimate the error y a − J a T I aq J q y qc in terms of a user-definable goal function Q : V a → R , (3.1)that is, we aim at estimating | Q ( y ac ) − Q ( J a T I aq J q y qc ) | . (3.2)We assume that Q is linear. Hence there exists some vector q ∈ V a such that Q ( y ) = q T y ∀ y ∈ V a . (3.3)We decompose the error (3.3) as follows: | Q ( y a ) − Q ( J a T I aq J q y qc ) | = | Q ( y a − y ac ) + Q ( y ac − J a T I aq J q y qc ) |≤ | Q ( y a − y ac ) | + | Q ( y ac − J a T I aq J q y qc ) | = | Q ( e a − ac ) | + | Q ( e ac − qc ) | (3.4)where e a − ac := y a − y ac , e ac − qc := y ac − J a T I aq J q y qc . (3.5)The first term | Q ( e a − ac ) | constitutes the modeling error and has been treated in [2]. The secondterm | Q ( e ac − qc ) | describes the error due to mesh coarsening and will be treated in the following.To facilitate the error analysis, we define the dual problems M ac g ac = q ,M qc g qc = J q T I aq T J a q . (3.6)We then have the basic dual identity for the goal-oriented error Q ( e ac − qc ) = q T ( y ac − J a T I aq J q y qc )= g ac T M ac ( y ac − J a T I aq J q y qc )= g ac T R ac ( J a T I aq J q y qc ) (3.7)with the primal residual given by R ac ( y ) := f ac − M ac y . (3.8) DAPTIVE REFINEMENT FOR THE QUASICONTINUUM APPROXIMATION 7
However, this quantity is too expensive to compute. We cannot solve for the dual solution g ac onthe atomistic scale, since this has the same computational complexity as solving for the uncoarsenedsolution y ac . To overcome this obstacle, we replace g ac by a dual solution from a coarser space.A first idea would be to use J a T I aq J q g qc instead of g ac . However, this turns out to be uselessdue to Galerkin orthogonality: Lemma 3.1 (Garlerkin Orthogonality) . We have J q T I aq T J a M ac ( y ac − J a T I aq J q y qc ) = 0 , (3.9) or equivalently J q T I aq T J a R ac ( J a T I aq J q y qc ) = 0 . (3.10) Proof.
Multiplying the equation for y ac from (2.24) by J q T I aq T J a from the left gives J q T I aq T J a M ac y ac = J q T I aq T J a f ac . (3.11)It is easy to see that J a J a T I aq J q = I aq J q . (3.12)Applying this to the equation for y qc from (2.24) leads to J q T I aq T J a M ac J a T I aq J q y qc = f qc . (3.13)Subtracting (3.13) from (3.11), substituting the definitions (2.26) of f ac and f qc and using (3.12)and (2.19) gives J q T I aq T J a M ac ( y ac − J a T I aq J q y qc )= J q T I aq T J a f ac − f qc = J q T I aq T J a J a T h − D T E ac D ( y bc − a a ) − K ( y bc − b a ) i − J q T I aq T h − D T E ac D ( I aq y bcq − a a ) − K ( I aq y bcq − b a ) i = 0 , (3.14)which completes the proof. (cid:3) Hence, replacing g ac in (3.7) by J a T I aq J q g qc always gives a zero estimate for the goal-orientederror. We need to use the dual solution from some space which is finer than V q to get a non-zeroestimate for the goal-oriented error, but which is coarser than V a to make it computable.To this end, we introduce an additional level of refinement and denote it as the partial continuum (pc) level. Similarly to the repatoms on the qc-level, we chose a set of 2 ¯ N pc-level repatoms whichis a subset of all atoms and a superset of the qc-repatoms. This means we choose indices ¯ ℓ ¯ for¯ = − ¯ N + 1 , . . . , ¯ N such that¯ ℓ − ¯ N +1 < ¯ ℓ − ¯ N +2 < · · · < ¯ ℓ ¯ < ¯ ℓ ¯ +1 < · · · < ¯ ℓ ¯ N − < ¯ ℓ ¯ N (3.15)and ¯ ℓ − ¯ N +1 = − M + 1 , ¯ ℓ − ¯ N +2 = − M + 2 , ¯ ℓ ¯ N − = M − , ¯ ℓ ¯ N = M. (3.16)To ensure that the pc-level is actually a refinement of the qc-level, we require every qc-level repatomto be a pc-level repatom. Hence there exist indices µ j such that ℓ j = ¯ ℓ µ j , j = − N + 1 , . . . , N. (3.17) DAPTIVE REFINEMENT FOR THE QUASICONTINUUM APPROXIMATION 8 ¯ ν − ¯ N +1 ¯ ν ¯ ¯ ν ¯ N − ν − N +1 ν j ν N − ac:pc:qc: − M +1 ℓ j = ¯ ℓ µ j M − ¯ N +1 ¯ = µ j ¯ N − N +1 j N Figure 3:
Levels of refinement: atomistic-continuum (ac), partial continuum (pc), and quasicontinuum (qc)with indices (above the respective chain) and intervals (below the respective chain).
Similar to the definition of ν j , we denote the number of atomistic intervals between two pc-levelrepatoms by ¯ ν ¯ := ¯ ℓ ¯ +1 − ¯ ℓ ¯ , ¯ = − ¯ N + 1 , . . . , ¯ N − . (3.18)See Figure 3 for an illustration of the three levels of refinement and the corresponding variables.We define the pc-level solution spaces V p := R N and V p := R N − , (3.19)the interpolation operators I ap ¯ ℓ ¯ +¯ k, ¯ := ¯ ν ¯ − ¯ k ¯ ν ¯ I ap ¯ ℓ ¯ +¯ k, ¯ +1 := ¯ k ¯ ν ¯ for ¯ = − ¯ N + 1 , . . . , ¯ N − , ¯ k = 0 , . . . , ¯ ν ¯ ,I api, ¯ := 0 otherwise ,I pqµ j + k,j := ν j − ¯ ℓ µ j + k + ¯ ℓ µ j ν j I pqµ j + k,j +1 := ¯ ℓ µ j + k − ¯ ℓ µ j ν j for j = − N + 1 , . . . , N − , k = 0 , . . . , µ j +1 − µ j ,I pq ¯ ,j := 0 otherwise , (3.20)the restriction operator R qpj, ¯ := δ µ j , ¯ for j = − N + 1 , . . . , N, ¯ = − ¯ N + 1 , . . . , ¯ N , (3.21)and the boundary operator J pi,j := δ ij for i = − ¯ N + 1 , . . . , ¯ N , j = − ¯ N + 3 , . . . , ¯ N − . (3.22)Note that the interpolation operators are defined in such a way that I aq factors as I aq = I ap I pq . (3.23) DAPTIVE REFINEMENT FOR THE QUASICONTINUUM APPROXIMATION 9 V q V p V a V q V p V a / / I pq o o R qp / / I ap * * I aq = I ap I pq O O J q O O J p O O J a (3.24) Figure 4:
Solution spaces and their operators for the ac-level, pc-level, and qc-level.
Altogether, all solution spaces and the corresponding operators are depicted in Figure 4.We solve the dual equation on the newly introduced partial continuum level: M pc g pc = J p T I ap T J a q , (3.25)where M pc := J p T I ap T ( D T E ac D + K ) I ap J p . (3.26)We use g pc to get an approximation of the goal-oriented error: Q ( y ac ) − Q ( J a T I aq J q y qc ) ≈ g pc T J p T I ap T J a R ac ( J a T I aq J q y qc ) . (3.27)Due to Galerkin Orthogonality (Lemma 3.1), we can subtract any vector in V q from the righthand side expression without changing its value. Thus, we subtract from g pc the linear interpolationof the nodal values of g pc at the qc-repatoms, J p T I pq R qp J p g pc , to get a vector in V p that is zeroat the qc-repatoms. This leads us to the error estimator η := (cid:2) g pc − J p T I pq R qp J p g pc (cid:3) T J p T I aq T J a R ac ( J a T I aq J q y qc )= (cid:2) g pc − J p T I pq R qp J p g pc (cid:3) T (cid:2) J p T I aq T J a f ac − M pc y pc (cid:3) (3.28)as an approximation to (3.7).In Section 5, we will construct a numerical method based on this error estimator. We will thencompute numerically how much the approximation (3.28) deviates from the precise error (3.7), andwe will determine how fine the partial refinement needs to be.4. Numerical Method
Based on the error estimator η , we need to decide what intervals at the qc-level shall be refined.To this end, we split η into individual contributions from each pc-level repatom: η = ¯ N − X ¯ = − ¯ N +3 η pc ¯ (4.1)where η pc ¯ := (cid:2) g pc − J p T I pq R qp J p g pc (cid:3) ¯ (cid:2) J p T I aq T J a R ac ( J a T I aq J q y qc ) (cid:3) ¯ . (4.2)Next, we map the values η pc ¯ to the qc-level intervals. Note that we already have η pc ¯ = 0 whenever¯ is a qc-repatom because we subtracted the qc-level interpolant from g pc in (3.28). We accumulate DAPTIVE REFINEMENT FOR THE QUASICONTINUUM APPROXIMATION 10 the remaining values η pc ¯ to the interval between the qc-repatoms j and j + 1 if the pc-repatom ¯ lies between the qc-repatoms j and j + 1, that is, µ j < ¯ < µ j +1 , to obtain η qcj := (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) µ j +1 − X ¯ = µ j +1 η pc ¯ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (4.3)Moreover, we have to decide how to compute the partial refinement. Here we subdivide each qc-interval ( ℓ j , ℓ j +1 ) into a fixed number Λ of subintervals of the same size, although other strategiesare possible, like letting Λ depend on the respective interval size ν j = ℓ j +1 − ℓ j , or dividing eachinterval into subintervals of different sizes. But here we stick to a fixed number Λ ∈ N and anequidistant subdivision. We employ the following algorithm for partial refinement:¯ ← − ¯ N + 1for j = − N + 1 , . . . , N − ω ← max(1 , ν j / Λ) σ ← σ ← σ < ν j ) σ ← min( σ + ω, ν j )¯ ν ¯ ← ⌊ σ − σ + ⌋ σ ← σ + ¯ ν ¯ ¯ ← ¯ + 1endendNote that this algorithm performs the necessary rounding if ν j is not divisible by Λ. If ν j ≤ Λ,then the interval gets fully refined up to the atomistic level.We then use this partial refinement algorithm and the global and local error estimators above toconstruct the following algorithm for adaptive mesh coarsening:(1) Start with the model fully coarsened in the continuum region.(2) Solve the primal problem (2.24) for y qc . (3) Determine the partial refinement pc according to the algorithm above.(4) Solve the dual problem (3.25) for g pc .(5) Compute the global error estimator η from (3.28).(6) If η ≤ τ gl , then stop.(7) Compute the local error estimators η qcj from (4.2) and (4.3).(8) Refine all intervals ( j, j + 1) with η qcj ≥ τ fac max k η qck (4.4)into two subintervals.(8) Go to (2).Here τ gl > τ fac = 10 is a reasonable choice.A possible improvement of step (8) would be to use the absolute value of the error estimator η qcj to decide into how many subintervals each interval ( j, j + 1) shall be refined, instead of justrefining into two subintervals. However, the magnitudes of the values η qcj differ considerably onlyin the first iterations. So this would only eliminate some of the early iterations which are still cheap DAPTIVE REFINEMENT FOR THE QUASICONTINUUM APPROXIMATION 11 iteration ν j max ν j η P N − j = − N +1 η qcj | Q ( e ac − qc ) | Table 1:
Advance of the algorithm until the error tolerance τ gl = 10 − is achieved. due to the small number of unknowns involved there, whereas the computationally expensive lateriterations will not be affected. Hence we expect the overall speedup to be small.5. Numerical Results
Now we present and discuss our numerical results. We choose boundary conditions y bcl = − M, y bcl = − M + 1 , y bcr = M − , y bcr = M. (5.1)The elastic moduli are given by k = 0 . k = 2, and k = 1, and the lattice constant is a = 1.We consider a chain of 4106 atoms, that is M = 2053. There is an atomistic region of 4 atomsfrom − − , − , , E a − ( I aq y qc ) and E a ( I aq y qc ) can be evaluated without interpolation, and we set and atoms − M +1 , − M +2 , M − , M to be continuum repatoms so that the boundary conditions can be set. Initially, the mesh in thecontinuum region is maximally coarsened, so that there are no more repatoms in addition to theones already mentioned. This gives N = 6 and 12 qc-level degrees of freedom. There are two largeelements of size ν − = ν = 2048, one on each side of the dislocation at the center of the chain,whereas the remaining elements necessarily have sizes ν j = 1 for j = − , − , − , − , , , , , y − y of the dislocation at the center of the chain,that is, the distance between the two atoms 0 and 1 to the left and right of the dislocation. Thecorresponding vector q ∈ V a reads as q = [0 , . . . , , − , , , . . . , T . (5.2)Table 1 shows how the algorithm proceeds for an error tolerance of τ gl = 10 − . One can easilysee how the elements get refined and the error drops. Also shown is the precise error | Q ( e ac − qc ) | which is available for this small model problem. Note that the column min ν j refers only to theelements in the continuum region which actually can be coarsened, excluding the padding aroundthe atomistic region and the boundary layer. DAPTIVE REFINEMENT FOR THE QUASICONTINUUM APPROXIMATION 12 −6 −4 −2 0 2 4 6124816326412825651210242048 m e s h s i z e −6 −4 −2 0 2 4 610 −24 −23 −22 −21 −20 −19 −18 −17 −16 −15 −14 −13 −12 −11 −10 −9 −8 −7 −6 −5 −4 −3 −2 η i q c −10 −5 0 5 10124816326412825651210242048 m e s h s i z e −10 −5 0 5 10 10 −24 −23 −22 −21 −20 −19 −18 −17 −16 −15 −14 −13 −12 −11 −10 −9 −8 −7 −6 −5 −4 −3 −2 η i q c −25 −20 −15 −10 −5 0 5 10 15 20 25124816326412825651210242048 m e s h s i z e −25 −20 −15 −10 −5 0 5 10 15 20 25 10 −24 −23 −22 −21 −20 −19 −18 −17 −16 −15 −14 −13 −12 −11 −10 −9 −8 −7 −6 −5 −4 −3 −2 η i q c Figure 5:
Mesh size ν j (green bar graph) and η qcj (black line graph) for error tolerances τ gl = 10 − (top), τ gl = 10 − (middle) and τ gl = 10 − (bottom). The bar graph in Figure 5 shows the size of the individual elements in the meshes before the firstiteration and after the error tolerances τ gl = 10 − and τ gl = 10 − are achieved. One can clearly seehow the elements get successively refined towards the center.The black lines in Figure 5 depict the decomposed error estimator η qcj . As a result of the Galerkinorthogonality, it vanishes whereever the mesh is fully refined up to the atomistic level. We can aswell read off the graphs how the algorithm tends to distribute the error uniformly over the wholemesh as the error tolerance is decreased. For τ gl = 10 − , we already achieved a close to flat regionfrom elements − . . . −
13 and 13 . . .
DAPTIVE REFINEMENT FOR THE QUASICONTINUUM APPROXIMATION 13 | Q ( e ac − qc ) | Λ η P η qcj η/ | Q ( e ac − qc ) | P η qcj / | Q ( e ac − qc ) | mesh 1 6.777614e-02 2 3.143618e-03 3.143618e-03 0.046382 0.0463824 8.351650e-03 8.351650e-03 0.123224 0.1232248 1.708501e-02 1.708501e-02 0.252080 0.252080 ∞ ∞ ∞ Table 2:
Efficiency of the error estimator, η . Table 2 shows the efficiency of the error estimator for different meshes and different values of Λ.Meshes 1, 2, and 3 refer to the meshes displayed in Figure 5. For comparison, we also include thevalue Λ = ∞ , which indicates that the pc-level mesh for the dual solution is fully refined to theatomistic level. One can read off from this table that the error indicator gets closer to the preciseerror if the mesh gets finer. For the coarse mesh 1, the actual error is considerably underestimated.However, this has only little impact on the final mesh since all elements have to be refined anyhow.For the finer meshes 2 and 3, η gets closer to the precise error | Q ( e ac − qc ) | .Also, we can see from Table 2 how the choice of Λ affects the error estimator. As expected, η gets closer to the precise error the larger Λ gets. For computational efficiency, we are interestedin keeping Λ small, though. We can read off from the table for the fine mesh 3 that already thesmallest possible value Λ = 2 gives a good estimate of the error, with a deviation of less than 20%.For Λ = 4, we already get a very precise estimate.The absolute accuracy of the error estimator is important, but even more important is how well η controls the mesh refinement, that is how efficient the resulting mesh is in terms of reaching aprescribed error tolerance with a minimal number of degrees of freedom. We now investigate theinfluence of the parameter Λ on the mesh quality.Table 3 shows the number of degrees of freedom and the corresponding error during the automaticmesh adaption for different values of Λ. The choices Λ = 4, Λ = 8, and Λ = ∞ lead to differentvalues of η , but all result in the same mesh. For this reason, these values share a common columnfor , , ∞ leads to slightly better mesh thatΛ = 2. However the difference is quite small so that the faster computation with Λ = 2 is very wellacceptable. References [1]
M. Ainsworth and J. T. Oden , A Posteriori Error Estimation in Finite Element Analysis , Wiley, New York,2000.[2]
M. Arndt and M. Luskin , Error estimation and atomistic-continuum adaptivity for the quasicontinuum ap-proximation of a Frenkel-Kontorova model , Multiscale Model. Simul., (2007), arXiv:0704.1924. Accepted.
DAPTIVE REFINEMENT FOR THE QUASICONTINUUM APPROXIMATION 14
Λ = 2 Λ = 4 Λ = 8 Λ = ∞ it | Q ( e ac − qc ) | η | Q ( e ac − qc ) | η η η Table 3:
Mesh efficiency. −10 −8 −6 −4 −2 p r e c i s e e rr o r Λ =2 Λ =4,8, ∞ Figure 6:
Mesh efficiency.
DAPTIVE REFINEMENT FOR THE QUASICONTINUUM APPROXIMATION 15 [3] ,
Goal-oriented atomistic-continuum adaptivity for the quasicontinuum approximation , Int. J. MultiscaleComput. Eng., (2007), arXiv:0708.0025. Accepted.[4]
S. Badia, M. L. Parks, P. B. Bochev, M. Gunzburger, and R. B. Lehoucq , On atomistic-to-continuum(AtC) coupling by blending , Technical report SAND 2007-5126J, Sandia National Laboratories, 2007.[5]
W. Bangerth and R. Rannacher , Adaptive Finite Element Methods for Differential Equations , Lectures inMathematics, ETH Z¨urich, Birkh¨auser, Basel, 2003.[6]
X. Blanc, C. Le Bris, and F. Legoll , Analysis of a prototypical multiscale method coupling atomistic andcontinuum mechanics , Math. Model. Numer. Anal., 39 (2005), pp. 797–826.[7]
W. A. Curtin and R. E. Miller , Atomistic/continuum coupling in computational materials science , ModellingSimul. Mater. Sci. Eng., 11 (2003), pp. R33–R68.[8]
M. Dobson and M. Luskin , Analysis of a force-based quasicontinuum approximation , Math. Model. Numer.Anal., (2007), arXiv:math.NA/0611543. To appear.[9]
W. E, J. Lu, and J. Z. Yang , Uniform accuracy of the quasicontinuum method , Phys. Rev. B, 74 (2006),p. 214115.[10]
W. E and P. Ming , Analysis of multiscale methods , J. Comput. Math., 22 (2004), pp. 210–219.[11]
J. Knap and M. Ortiz , An analysis of the quasicontinuum method , J. Mech. Phys. Solids, 49 (2001), pp. 1899–1923.[12]
P. Lin , Theoretical and numerical analysis for the quasi-continuum approximation of a material particle model ,Math. Comput., 72 (2003), pp. 657–675.[13] ,
Convergence analysis of a quasi-continuum approximation for a two-dimensional material without defects ,SIAM J. Numer. Anal., 45 (2007), pp. 313–332.[14]
W. K. Liu, E. G. Karpov, S. Zhang, and H. S. Park , An introduction to computational nanomechanics andmaterials , Comput. Methods Appl. Mech. Eng., 193 (2004), pp. 1529–1578.[15]
M. Marder , Condensed Matter Physics , John Wiley & Sons, 2000.[16]
R. Miller and E. Tadmor , The quasicontinuum method: Overview, applications and current directions , J.Comput. Aided Mater. Des., 9 (2002), pp. 203–239.[17]
R. Miller, E. B. Tadmor, R. Phillips, and M. Ortiz , Quasicontinuum simulation of fracture at the atomicscale , Modelling Simul. Mater. Sci. Eng., 6 (1998), pp. 607–638.[18]
J. T. Oden, S. Prudhomme, and P. Bauman , On the extension of goal-oriented error estimation and hierar-chical modeling to discrete lattice models , Comput. Methods Appl. Mech. Eng., 194 (2005), pp. 3668–3688.[19] ,
Error control for molecular statics problems , Int. J. Multiscale Comput. Eng., 4 (2006), pp. 647–662.[20]
C. Ortner and E. S¨uli , A-posteriori analysis and adaptive algorithms for the quasicontinuum method in onedimension , Research Report NA-06/13, Oxford University Computing Laboratory, 2006.[21] ,
A-priori analysis of the quasicontinuum method in one dimension , Research Report NA-06/12, OxfordUniversity Computing Laboratory, 2006.[22]
M. L. Parks, P. B. Bochev, and R. B. Lehoucq , Connecting atomistic-to-continuum coupling and domaindecomposition , Technical report SAND 2007-0704J, Sandia National Laboratories, 2007.[23]
V. B. Shenoy, R. Miller, E. B. Tadmor, D. Rodney, R. Phillips, and M. Ortiz , An adaptive finiteelement approach to atomic-scale mechanics - the quasicontinuum method , J. Mech. Phys. Solids, 47 (1999),pp. 611–642.[24]
E. B. Tadmor, M. Ortiz, and R. Phillips , Quasicontinuum analysis of defects in solids , Philos. Mag. A, 73(1996), pp. 1529–1563.[25]