Numerical scheme for solving the nonuniformly forced cubic and quintic Swift-Hohenberg equations strictly respecting the Lyapunov functional
NNumerical scheme for solving the nonuniformly forced cubic and quinticSwift-Hohenberg equations strictly respecting the Lyapunov functional
D. L. Coelho a, ∗ , E. Vitral b , J. Pontes a , N. Mangiavacchi a a GESAR Group/UERJ - State University of Rio de Janeiro, 20940-903 Rua Fonseca Teles 121, Rio de Janeiro, RJ, Brazil b Aerospace Engineering and Mechanics Department, University of Minnesota, 107 Akerman Hall, Minneapolis, MN, USA55455-0153
Abstract
Computational modeling of pattern formation in nonequilibrium systems is a fundamental tool for studyingcomplex phenomena in biology, chemistry, materials science and engineering. The pursuit for theoreticaldescriptions of some among those physical problems led to the Swift-Hohenberg equation (SH3) which de-scribes pattern selection in the vicinity of instabilities. A finite differences scheme, known as StabilizingCorrection (Christov & Pontes; 2001 DOI: 10.1016/S0895-7177(01)00151-0), developed to integrate the cubicSwift-Hohenberg equation in two dimensions, is reviewed and extended in the present paper. The originalscheme features Generalized Dirichlet boundary conditions (GDBC), forcings with a spatial ramp of the con-trol parameter, strict implementation of the associated Lyapunov functional, and second order representationof all derivatives. We now extend these results by including periodic boundary conditions (PBC), forcingswith gaussian distributions of the control parameter and the quintic Swift-Hohenberg (SH35) model. Thepresent scheme also features a strict implementation of the functional for all test cases. A code verificationwas accomplished, showing unconditional stability, along with second order accuracy in both time and space.Test cases confirmed the monotonic decay of the Lyapunov functional and all numerical experiments exhibitthe main physical features: highly nonlinear behaviour, wavelength filter and competition between bulk andboundary effects.
Keywords:
Pattern Formation, Nonlinear Systems, Swift-Hohenberg Equation, Finite Difference Methods
1. Introduction
Understanding the selection and orientation mechanisms of the spatial structures, their symmetries andinstabilities represent a major theme of theoretical and experimental research in computational modeling ofself-organization phenomena. Nowadays, spatiotemporal organization is recognized to be related to severaltechnological problems in physics, chemistry [1], nonlinear optics [2], and materials science [3, 4, 5]. Forcomputational material science and engineering a long-term goal is to apply such modeling and numericalframework to multiple scales and increasingly larger scale systems in order to achieve predictive capabilities [6, ∗ Corresponding author. Tel.: +55 (21) 2332-4733
Email addresses: [email protected] (D. L. Coelho), [email protected] (D. L. Coelho)
August 3, 2020 a r X i v : . [ n li n . PS ] J u l , 8, 9].The pursuit for theoretical descriptions of some among those physical problems led to the Swift-Hohenberg(SH) equation [10]. It is a widely accepted model for describing pattern formation in many physical systemspresenting symmetry breaking instabilities, such as the emergence of convection rolls in a thin layer of fluidheated from below. It first appeared in the framework of Bénard thermal convection between two “infinite”horizontal surfaces with distinct temperatures, where J. Swift and P. Hohenberg [10] performed a reductionof the full dynamics, led by the Boussinesq equations, to the slow modes dynamics. Such reduction wasalso accomplished for reaction-diffusion systems where a similar model could be derived with an additionalquadratic nonlinearity (SH23) [3, 11]. More recently, an extension of the original Swift-Hohenberg equation(SH3) was proposed with the introduction of a destabilizing cubic term and a quintic one. The resultingquintic SH equation (SH35) admits the coexistence of stable uniform and structured solutions, which maycoexist and extends the use of the Swift-Hohenberg equation to the study of localized structures [12].Roughly speaking, the SH equation has three basic pattern-forming mechanisms, related to each term ofthe equation. The linear term is accompanied by a control parameter (related to temperature in the caseof Rayleigh Bénard convection) and will be addressed as a forcing parameter. It represents how far fromthe onset of the instability the system is. Secondly, the wavelength filter contains all the equation spatialderivatives and tends to lead the pattern to have the chosen critical wavelength. Finally the nonlinear termsplay the interaction between the modes involved in the dynamics. In the SH3, the cubic term is responsiblefor saturating the linear growth. In the SH35, the cubic term is destabilize the dynamics while the quinticone is responsible for the saturation. Some theoretical works explored nonuniform forcings for the case ofRayleigh Bénard convection [13, 14].The SH model has a potential form such that we can associate its dynamics to a free energy functional,usually called the Lyapunov functional. In the last years, this potential nature has been explored in thematerials science framework as a physical feature associated with the system’s free energy and can be relatedto the description of materials patterning. In this context, the SH model is part of the phase-field theory,which was developed from statistical mechanics principles, and whose goal is to obtain governing equationsfor microstructure evolution; it connects thus thermodynamic and kinetic properties with microstructure viaa mathematical formalism [5, 7, 15, 16]. The phase-field theory is a descendant of the van der Waals, Cahn-Hilliard and Landau type classical field theoretical approaches, originating from a single order parametergradient theory of Langer [17, 18]. In the phase-field theory, the local state of matter is characterized byan order parameter ψ ( x , t ) , called the phase-field variable, which monitors the transition between phases ofdistinct order of property. It can represent the structural order parameter that measures local crystallinity,the composition of a phase, the degree of a molecular ordering, etc. In the case of the phase-field crystal(PFC) approach, a conserved dynamics is assumed for this order parameter. Although this formalism willnot be addressed in this paper, it is worth mentioning that many interesting numerical works explored theSH equation in a PFC framework [19, 20, 21]. 2here are many works presenting fruitful numerical frameworks for the Swift-Hohenberg equation. Inthe beginning of this endeavor, Greenside & Coughran Jr. (1984)[22] achieved relevant results through a fi-nite differences approach by using a backward Euler time-integration scheme with rigid [23, 24] and periodicboundary conditions (PBC). Those types of boundary conditions were also explored by Manneville (1990)[25]and Cross (1994)[26]. These works considered uniform control parameters. Other authors employed nonuni-form distributions (ramps and gaussians) for the control parameter and presented very interesting results aswell, displaying rich competition between bulk and boundary effects. A ramped control parameter mostlyappears in the context of rigid boundaries [27, 28], which will be addressed here as Generalized Dirichletboundary conditions (GDBC), following [27, 29, 28].Spectral methods have become largely used due to their highly spatial accurate solutions [30, 31, 32].They usually treat the nonlinear terms explicitly in the physical space and therefore are called pseudo-spectral methods. Fully explicit schemes leads to severe stability restrictions on the time step, because of thefourth-order nonlinear partial differential equation originated from the non-conserved dynamics. In the caseof the phase-field crystal (PFC) approach (conserved dynamics), one derives a six-order nonlinear partialdifferential equation, which is even more restrictive for the time step selection. While PFC is not a subjectof this work, our proposed scheme can be readily translated into a numerical framework for many PFCproblems.It is also pertinent to point out that operator splitting methods enable parallel code processing imple-mentation, which allows faster simulations. The usage of sparse matrices also plays a huge role in loweringcomputational cost in the simulations, and therefore it was implemented in the present work.The present article is a continuation of the work presented by Christov & Pontes (2002)[29] where thoseauthors employed a finite differences scheme originally introduced by Douglas & Rachford in the frameworkof the temperature equation, for solving the cubic Swift-Hohenberg equation (SH3) in two dimensions, withuniform forcing, GDBC, and strict implementation of the associated functional. The scheme is known asStabilizing Correction [33, 34]. Here, we expand this study for additional nonlinearities, boundary conditionsand a spatially dependent control parameter. This includes the quintic version of the Swift-Hohenbergequation (SH35) with both GDBC and PBC, nonuniform gaussian distributions of the control parameterin addition to ramps, and also a strict implementation of the associated Lyapunov functional. The schemefeatures a semi-implicit time splitting with second order representation of time and space derivatives. Someother works also used stabilizing correction scheme for similar fourth-order nonlinear differential equations,such as the anisotropic damped Kuramoto-Sivashinsky equation [35].A code verification is performed to confirm the unconditional stability of the scheme, along with its secondorder convergence in both time and space. Since nonlinear partial differential equations do not easily haveanalytical solutions, which is the SH equation case, the adopted method to verify the order of accuracy ofthe code is the method of manufactured solutions (MMS). It provides a convenient way of verifying theimplementation of nonlinear numerical algorithms, since we can use an artificial (manufactured) solution for3his purpose [36, 37, 35].This paper is organized as follows. Section 2 briefly discusses the properties of the SH equation andthe numerical scheme for solving the cubic and quintic versions of the equation with GDBC and PBC,nonuniform forcings, and the discrete implementation of the associated Lyapunov functional. Section 3contains all verification procedures adopted. Section 4 presents the numerical experiments. The results aregrouped in three parts; the first one showing the steady state patterns obtained in ten simulations withthe SH3 equation; the second one showing six simulations with the SH35 and the last one with discussions.Section 5 contains the final remarks of the work.
2. Mathematical modelling and numerical scheme
The mathematical description of the dimensionless governing equations is briefly discussed in the beginningof this section. Then, the details of the numerical framework are exposed as an extension of Christov &Pontes (2002)[29] work. For both SH3 and SH35 discretizations, a finite difference semi-implicit scheme ofsecond order accuracy in time and space is presented. Finally, the choices for the operator splitting method,spatial discretizations and mesh types are clarified.
The SH equation has the so-called gradient dynamics, which means there is a potential function of theorder parameter ψ ( x , t ) , (known as a Lyapunov functional), that has the property of decreasing monotonicallyduring the evolution [29]. This dynamic equation can be derived from the L -gradient flow of the Lyapunovenergy functional: F [ ψ ] := (cid:90) Ω d x (cid:26) − (cid:15) ( x ) ψ + α (cid:2) ( q + ∇ ) ψ (cid:3) − β ψ + γ ψ (cid:27) ; (1) ∂ F [ ψ ] ∂t = − (cid:90) Ω d x (cid:18) ∂ψ∂t (cid:19) ≤ . (2)where Ω represents the domain whose size is commensurate with the length scales of the patterns. Weconsider a regular domain Ω : { x ∈ [0 , L x ] , y ∈ [0 , L y ] } . Also, Eq. 2 denotes the nonincreasing behaviorof the Lyapunov functional, which monotonically decreases until steady state is reached [29]. By taking thevariational derivative of Eq. 1 in L norm, the following Swift-Hohenberg equation is obtained: ∂ψ∂t = − δ F [ ψ ] δψ = (cid:15) ( x ) ψ − α (cid:0) q + ∇ (cid:1) ψ + βψ − γψ . (3)As discussed by Christov et al. [29, 27], the appropriate boundary conditions for this system shouldconsider only production and dissipation in the bulk and not on the boundary ( ∂ Ω ). This can be achievedby assuming Generalized Dirichlet boundary conditions (GDBC), v = ∂v/∂n = 0 , where n stands for theoutward normal direction to the boundary. Periodic boundary conditions (PBC) will also be addressed inthis work, since the absence of boundaries also leads boundary ( ∂ Ω ) integrals to vanish.4e expand Eq. 3 by considering the following laplacian operator in cartesian coordinates: ∇ ≡ ∂ /∂x + ∂ /∂y , we have: ∂ψ∂t = (cid:15) ( x ) ψ − αq ψ − αq ∇ ψ − α ∇ ψ + βψ − γψ = (cid:15) ( x ) ψ − αq ψ − αq ∂ ψ∂x − αq ∂ ψ∂y − α ∂ ψ∂x − α ∂ ψ∂x ∂y − α ∂ ψ∂y + βψ − γψ . (4) Table 1: Parameters assumed for the governing equations studied in this work.
Equation Nonlinearity α β γ q SH3 cubic 1.0 -1.0 0.0 1.0SH35 quintic 1.0 1.0 1.0 1.0By using the parameters and values of Tab. 1 in Eq. 3, we define the forms of the SH3 and SH35. Thediscretization of these governing equations is described in the following subsection.
In order to construct a Crank-Nicolson second order in time numerical scheme, we adopt the followingrepresentation proposed by Christov and Pontes (2002) [29] for the time derivative of Eq. 3, where the RHSis evaluated at the middle of the time step ∆ t . The updated scheme, now including a quintic term takes theform: ψ n +1 − ψ n ∆ t = (cid:20) (cid:15) ( x ) − αq − αq ∂ ∂x − αq ∂ ∂y − α ∂ ∂x − α ∂ ∂x ∂y − α ∂ ∂y + β (cid:0) ψ n +1 (cid:1) + ( ψ n ) − γ (cid:0) ψ n +1 (cid:1) + ( ψ n ) (cid:35) (cid:18) ψ n +1 + ψ n (cid:19) . (5)The superscript ( n + 1) refers to the next time to be evaluated and n , to the current one. The parametervalues for α , β , γ and q are chosen according to Tab. 1.The RHS terms of Eq. 5 are grouped in three parts, the first and the second ones containing the semi-implicit operators Λ n +1 / x and Λ n +1 / y that will act on the variable (cid:0) ψ n +1 + ψ n (cid:1) / , and a function f n +1 / ,evaluated at the middle of the time step ∆ t . This function will contain explicit terms only, in the final discreteform of Eq. 5. Space derivatives are represented by centered second order formulæ. Implicit terms are chosenhaving in mind to construct negative definite operators that will assure the stability of the scheme. Thefactor / multiplying (cid:0) ψ n +1 + ψ n (cid:1) in the RHS of the above equation is included in the operators Λ n +1 / x and Λ n +1 / y , leading to the following target scheme: ψ n +1 − ψ n ∆ t = (cid:16) Λ n +1 / x + Λ n +1 / y (cid:17) (cid:0) ψ n +1 + ψ n (cid:1) + f n +1 / . (6)5or the SH3 equation, the operators Λ n +1 / x , Λ n +1 / y and f n +1 / are defined as: Λ n +1 / x = 12 (cid:34) − α (cid:18) ∂ ∂x + q (cid:19) − β (cid:0) ψ n +1 (cid:1) + ( ψ n ) (cid:35) ;Λ n +1 / y = 12 (cid:34) − α (cid:18) ∂ ∂y + q (cid:19) − β (cid:0) ψ n +1 (cid:1) + ( ψ n ) (cid:35) ; f n +1 / = 12 (cid:20) (cid:15) ( x ) − α (cid:18) q ∂ ∂x + 2 q ∂ ∂y + 2 ∂ ∂x ∂y (cid:19)(cid:21) (cid:0) ψ n +1 + ψ n (cid:1) , (7)and for the SH35: Λ n +1 / x = 12 (cid:34) − α (cid:18) ∂ ∂x + q (cid:19) − γ (cid:0) ψ n +1 (cid:1) + ( ψ n ) (cid:35) ;Λ n +1 / y = 12 (cid:34) − α (cid:18) ∂ ∂y + q (cid:19) − γ (cid:0) ψ n +1 (cid:1) + ( ψ n ) (cid:35) ; f n +1 / = 12 (cid:20) (cid:15) ( x ) − α (cid:18) q ∂ ∂x + 2 q ∂ ∂y + 2 ∂ ∂x ∂y (cid:19) + β (cid:0) ψ n +1 (cid:1) + ( ψ n ) (cid:35) (cid:0) ψ n +1 + ψ n (cid:1) . (8) Since the operators Λ n +1 / x , Λ n +1 / y and the function f n +1 / in Eqs. 6, 7 and 8 contain implicit termsdue to the time discretization and linearization of the nonlinear terms, we do internal iterations. They arerequired to secure the approximation of the nonlinearities in the scheme (Eq. 6) at each time step. Theiterations loop proceeds until convergence is attained by monitoring the L ∞ norm. The internal iterationsscheme reads: ψ n +1 ,p +1 − ψ n ∆ t = (cid:16) Λ n +1 / ,px + Λ n +1 / ,py (cid:17) (cid:0) ψ n +1 ,p +1 + ψ n (cid:1) + f n +1 / ,p , (9)where the index ( p ) refers to the internal iteration number. The superscript ( n + 1 , p + 1) identifies the newiteration, while ( n ) are the values of the previous time step. The superscript ( n + 1) for the nonlinear termin the function f n +1 / will be replaced by ( n, p ) , which stands for the values obtained from the previousiteration.The operators Λ n +1 / ,px , Λ n +1 / ,py function f n +1 / ,p are redefined as follows, for the SH3 equation: Λ n +1 / ,px = 12 (cid:34) − α (cid:18) ∂ ∂x + q (cid:19) − β (cid:0) ψ n +1 ,p (cid:1) + ( ψ n ) (cid:35) , Λ n +1 / ,py = 12 (cid:34) − α (cid:18) ∂ ∂y + q (cid:19) − β (cid:0) ψ n +1 ,p (cid:1) + ( ψ n ) (cid:35) ,f n +1 / ,p = 12 (cid:20) (cid:15) ( x ) − α (cid:18) q ∂ ∂x + 2 q ∂ ∂y + 2 ∂ ∂x ∂y (cid:19)(cid:21) (cid:0) ψ n +1 ,p + ψ n (cid:1) , (10)and for the SH35: 6 n +1 / ,px = 12 (cid:34) − α (cid:18) ∂ ∂x + q (cid:19) − γ (cid:0) ψ n +1 ,p (cid:1) + ( ψ n ) (cid:35) , Λ n +1 / ,py = 12 (cid:34) − α (cid:18) ∂ ∂y + q (cid:19) − γ (cid:0) ψ n +1 ,p (cid:1) + ( ψ n ) (cid:35) ,f n +1 / ,p = 12 (cid:20) (cid:15) ( x ) − α (cid:18) q ∂ ∂x + 2 q ∂ ∂y + 2 ∂ ∂x ∂y (cid:19) + β (cid:0) ψ n +1 ,p (cid:1) + ( ψ n ) (cid:35) (cid:0) ψ n +1 ,p + ψ n (cid:1) . (11)The iterations proceed until the following criterion for the L ∞ norm is satisfied with δ = 1 . × − : L ∞ = max | ψ n +1 ,p +1 − ψ n +1 ,p | max | ψ n +1 ,p +1 | ≤ δ ; (12)so that the last iteration gives the sought function ψ in the new time ψ n +1 ≡ ψ n +1 ,p +1 . Although employing sparse matrices for the operators, computationally solving Eq. 9 still represents acostly procedure. In order to reduce such computational effort and errors (from discretization and floating-point operations), the operators splitting method was adopted. The splitting of Eq. 6 is made according tothe Douglas second scheme (also known as scheme of stabilizing correction, shown by [27],[29]), and is brieflyreviewed. The following equations represent a consistent approximation of the original scheme: (cid:101) ψ − ψ n ∆ t = Λ n +1 / ,px (cid:101) ψ + Λ n +1 / ,py ψ n + f n +1 / ,p + (Λ n +1 / ,px + Λ n +1 / ,py ) ψ n , (13) ψ n +1 ,p +1 − (cid:101) ψ ∆ t = Λ n +1 / ,py ( ψ n +1 ,p +1 − ψ n ) , (14)where (cid:101) ψ is an intermediary estimation of ψ at the new time step. In order to show that the splitting representsthe original scheme, we rewrite Eqs. 13 and 14 in the form: (cid:16) E − ∆ t Λ n +1 / ,px (cid:17) (cid:101) ψ = (cid:16) E + ∆ t Λ n +1 / ,px (cid:17) ψ n + 2∆ t Λ n +1 / ,py ψ n + ∆ tf n +1 / ,p , (15) (cid:16) E − ∆ t Λ n +1 / ,py (cid:17) ψ n +1 ,p +1 = (cid:101) ψ − ∆ t Λ n +1 / ,py ψ n , (16)where E is the unity operator. Rearranging these equations, the intermediate variable (cid:101) ψ is eliminated andthe result may be rewritten as: (cid:16) E − ∆ t Λ n +1 / ,px (cid:17) (cid:16) E − ∆ t Λ n +1 / ,py (cid:17) ψ n +1 ,p +1 = (cid:16) E + ∆ t Λ n +1 / ,px (cid:17) ψ n + 2∆ t Λ n +1 / ,py ψ n ++ ∆ tf n +1 / ,p − (cid:16) E − ∆ t Λ n +1 / ,px (cid:17) ∆ t Λ n +1 / ,py ψ n . (17)This result may be rewritten as: 7 E + ∆ t Λ n +1 / ,px Λ n +1 / ,py (cid:17) ψ n +1 ,p +1 − ψ n ∆ t = (Λ n +1 / ,px + Λ n +1 / ,py )( ψ n +1 ,p +1 + ψ n ) + f n +1 / ,p , (18)where E is the unity operator. A comparison with Eq. 6 shows that Eq. 18 is actually equivalent to theformer except by the positive definite operator having a norm greater than one: B ≡ E + ∆ t Λ n +1 / ,px Λ n +1 / ,py = E + O (∆ t ) , (19)which acts on the term ( ψ n +1 ,p +1 − ψ n ) / ∆ t . This means that this operator does not change the steady statesolution. Furthermore, since || B || > , the scheme given by Eqs. 15, and 16 is more stable than the targetone (Eq. 6). We solve numerically the SH equation with GDBC and PBC. For the case of rigid boundary conditionswe adopt a staggered grid, as illustrated in Fig. 1. Consider a staggered mesh in both spatial directions,namely x i = − ∆ x i ∆ x, ∆ x ≡ L x n x − , y j = − ∆ y j ∆ y, ∆ y ≡ L y n y − , where n x and n y are the number of points in x -and y -directions, respectively. The mesh pattern is shown inFig. 1. Let ψ i,j be an arbitrary set function defined on the above described mesh. We confine ourselves to thecase of constant coefficients. The PBC, which is less restrictive than the GDBC, borrows the crystallographyconcept of unit cell, which is a repeating pattern representative of the material. That is, using PBC we havea domain that would also act as a unit cell for a infinite two-dimensional surface, with cells being periodicallyrepeated all around the bidimensional domain. ψ ( x , ) = ∂ψ∂ y (cid:12)(cid:12)(cid:12)(cid:12) y = = ψ ( , y ) = ∂ψ∂ y (cid:12)(cid:12)(cid:12)(cid:12) x = = ψ ( L x , y ) = ∂ψ∂ y (cid:12)(cid:12)(cid:12)(cid:12) x = L x = ψ ( x , L y ) = ∂ψ∂ y (cid:12)(cid:12)(cid:12)(cid:12) y = L y = ◦ ◦◦ ◦ ◦◦◦ ◦◦ ◦◦◦ ◦◦◦◦ j = n y j = n y − j = j = i = n x i = n x − i = i = ψ nx , ψ , ψ , ··· ψ nx − , ψ nx , ψ , ψ nx , ny ψ , ny ψ , ny ··· ψ nx − , ny ψ nx , ny ψ , ny ψ nx , ny − ψ , ny − ψ , ny − ··· ψ nx − , ny − ψ nx , ny − ψ , ny − ... ... ... ... ... ... ... ψ nx , ψ , ψ , ··· ψ nx − , ψ nx , ψ , ψ nx , ψ , ψ , ··· ψ nx − , ψ nx , ψ , ψ nx , ny ψ , ny ψ , ny ··· ψ nx − , ny ψ nx , ny ψ , ny (a) (b) Figure 1: The grids and boundary conditions used in this work. (a) The “
Staggered ” grid. In the case of GDBC, values of ψ inthe first two lines and in the two last ones, and also in the first two columns and in the last two ones are assigned to zero. (b)The periodic domain. ψ at the points of a uniform grid. Considering thefunction ψ in a bidimensional domain, we can define the derivatives by truncating the Taylor expansion andarranging the equations.Second derivatives with second order accuracy can be written as: ∂ ψ i,j ∂x = ψ i,j − − ψ i,j + ψ i,j +1 ∆ x ; ∂ ψ i,j ∂y = ψ i − ,j − ψ i,j + ψ i +1 ,j ∆ y . (20)And fourth derivatives with second order accuracy are written as: ∂ ψ i,j ∂x = ψ i − ,j − ψ i − ,j + 6 ψ i,j − ψ i +1 ,j + ψ i +2 ,j ∆ x ,∂ ψ i,j ∂y = ψ i,j − − ψ i,j − + 6 ψ i,j − ψ i,j +1 + ψ i,j +2 ∆ y ,∂ ψ i,j ∂x ∂y = 1∆ x ∆ y (cid:32) ψ i − ,j − − ψ i,j − + ψ i +1 ,j − − ψ i − ,j + 4 ψ i,j − ψ i +1 ,j ++ ψ i − ,j +1 − ψ i,j +1 + ψ i +1 ,j +1 (cid:33) . Standard second order representations of spatial derivatives are adopted in uniform and structured grids.In order to verify the correctness of the implementation, we conducted convergence tests using the methodof manufactured solutions (MMS). The results can be found in Sec. 3.
3. Code verification
One important issue concerning the simulations is the time and mesh size selection, such that we seekreasonable choices inside the stable region of the proposed numerical scheme. The governing equation 3 canbe rewritten in the form
P ψ ( x , t ) = 0 , where P is an operator containing all the partial derivatives and termsacting on the order parameter ψ ( x , t ) . The consistency of the numerical scheme can be easily verified since P ψ ( x , t ) − P ∆ t, ∆ x, ∆ y ψ ( x , t ) −→ as ∆ t, ∆ x, ∆ y −→ , where P ∆ t, ∆ x, ∆ y is the finite difference discretizationof P . In this section, scheme stability, free energy functional decay and convergence tests were performed toverify the code implementation. Here, following Christov and Pontes (2001)[29], we study the effect of the time step in the structureevolution by assessing the rate of change in time of the pattern during the simulation. We do this bymonitoring ˙ L , the relative norm rate of change defined as: ˙ L = 1∆ t n x (cid:80) i =1 n y (cid:80) j =1 | ψ n +1 i,j − ψ ni,j | n x (cid:80) i =1 n y (cid:80) j =1 | ψ n +1 i,j | , (21)9hich roughly corresponds to the ratio between the spatial average of the modulus of time derivative ∂ψ/∂t and the spatial average of the modulus of the function itself. The calculations begin from a random initialcondition and proceeded until ˙ L ≤ × − , which is our criterion for reaching the steady state. Followingthis implementation, Fig. 2 shows the system state at t = 100 and the steady state attained in six simulationsrun with different time steps ∆ t . Fig. 3 shows the evolution of the associated ˙ L and the curves of theaccomplished internal iterations at each time step. This group of simulations was run with GDBC. t = 100 t = 100 t = 100 t = 100 t = 100 t = 100∆ t = 0 .
01 ∆ t = 0 . t = 0 . t = 1 . t = 2 . t = 10 . t = 3062 t = 3066 t = 3142 t = 3412 t = 4522 t = 30230∆ t = 0 .
01 ∆ t = 0 . t = 0 . t = 1 . t = 2 . t = 10 . Figure 2: Pattern developed with the SH3 model, GDBC, forced with a spatial ramp of the control parameter given by . ≤ (cid:15) ( x ) ≤ . , and six different time steps. Top row: transient states at t = 100 . Bottom: the steady state for the six timesteps. All tests started from the exactly same initial condition (pseudo-randomly generated) for a × nodes domain. Notethat the use of larger time steps results in smaller number of required steps to attain the same “time” ( t = 100 for instance).Errors with larger time steps result in delay in the emergence of the pattern and thus later steady states is expected. ˙ L [-] i ≈ ii iii iv vvii. ∆ t = ∆ t = ∆ t = ∆ t = ∆ t = ∆ t = I n t e r n a lit e r a ti on s (-) iiiiii iv vvii. ∆ t = ∆ t = ∆ t = ∆ t = ∆ t = ∆ t = ˙ L [-] i ≈ ii iii iv vvii. ∆ t = ∆ t = ∆ t = ∆ t = ∆ t = ∆ t = I n t e r n a lit e r a ti on s (-) iiiiii iv vvii. ∆ t = ∆ t = ∆ t = ∆ t = ∆ t = ∆ t = (a) (b) Figure 3: Scheme stability analysis from the numerical integration of the Swift-Hohenberg equation. (a) ˙ L and (b) internaliterations are shown as a function of time t for a × nodes domain for all six tests are presented in Fig. 2. The same steadypattern is reached for all six analyzed cases independently of the time step. The simulations started from the same randominitial conditions. ∆ t = 0 . , which presents the same ˙ L decaycurve as for ∆ t = 0 . . The choice for the remainder of this work was ∆ t = 0 . for the SH3 simulationsand ∆ t = 0 . for the SH35, so that lower computational time is required without losing physically consistenttransient results given the numerical scheme stability. Table 2: Computational time spent on each of the simulations presented in 2. simulations.
Time Minimum Maximum Steady Computationalstep iterations iterations state time spent ∆ t = 0 . t = 3062
03 weeks ∆ t = 0 . t = 3066 ∆ t = 0 . t = 3142 ∆ t = 1 . t = 3412 ∆ t = 2 . t = 4522 ∆ t = 10 . t = 30230 The adopted method to verify the order of accuracy of the code is the method of manufactured solutions(MMS). It provides a convenient way of verifying the implementation of nonlinear numerical algorithms byusing a manufactured (artificial) solution for such purpose [36, 37, 35]. In terms of the proposed problem, wetake all the members of the SH equation 3 and consider the following differential equation: F ( ψ ) = 0 , (22)where ψ is the order parameter function that satisfies Eq. 22 and therefore is the PDE solution. TheMMS consists of adopting an arbitrary function to be the manufactured solution, ψ m ( x , t ) , and since thisfunction is not likely to solve the PDE, a source term is expected, s m . This term can be seen as an additionalforcing function, leading to a modified operator with this new source: ¯ F ( ψ ) ≡ F ( ψ ) − s m . (23)For the previous equation, ¯ F ( ψ m ) = 0 and ¯ F ( ψ ) = − s m . Following this new approach to the problem, wefind an approximate numerical solution, ψ k , for the discretized problem so that ¯ F ( ψ k ) = 0 or F ( ψ k ) = s m .This source term is a minimal intrusion to the code’s formulation. The chosen function is periodic with awavenumber q and is defined as: ψ m ( x , t ) = ψ + ψ xy cos[ q ( x + y )] e at , (24)11here all parameters employed in the manufactured solution and in the differential equation are present onTab. 3. Table 3: Parameters assumed for the manufactured solution adopted in the convergence analysis.
MS Parameters Value SH Parameters Value q . √ q (cid:15) -0.1 q ∗ √ q q ψ α ψ xy (cid:112) | (cid:15) | β -1.0 a γ L norm, defined as follows: L = n x (cid:80) i =1 n y (cid:80) j =1 | ( ψ k ) ni,j − ( ψ m ) ni,j | n x (cid:80) i =1 n y (cid:80) j =1 | ( ψ m ) ni,j | / . (25)In the previous section, a second-order scheme was presented and therefore the formal order of accuracy istwo. The observed order of accuracy of the code can be acquired from the global discretization error formeshes with different grid spacing, and can described by the following relation: p = ln (cid:0) L B /L A (cid:1) ln( r ) , (26)where L A and L B are the L norm for meshes A and B respectively, and r is the ratio of the grid resolution, g r , given in number of mesh nodes per critical wavelength of mesh B to A. The meshes employed for thepresent tests, the number of nodes, and grid resolutions are shown in Tab. 4, and the L curves for thenumerical experiments can be seen in Fig. 4. 12 able 4: Meshes employed for the code verification by MMS. The chosen domain for the test has × critical wavelengths. Theobserved order of accuracy of the code in Fig. 4 can be calculated using Eq. (26) for two cases: q = 0 . √ q and q ∗ = √ q .In the first case the formal order of accuracy is achieved, p ≈ . , while in the second case p ∗ ≈ . due to discretization errorcanceling. Mesh Number Grid Convergence rate ( p ) of nodes resolution ( g r ) (cid:0) q = 0 . √ q (cid:1) (cid:0) q ∗ = √ q (cid:1) Mesh A × × . . Mesh C × . . Mesh D × . . Mesh E × . . g r [-]1e-011e-021e-031e-041e-051e-061e-071e-08 L [-] L norm2 nd order slope4 th order slope L [-] g r = g r = g r = g r = g r = (a) (b) g r [-]1e-011e-021e-031e-041e-051e-061e-071e-08 L [-] L norm2 nd order slope4 th order slope L [-] g r = g r = g r = g r = g r = (c) (d) Figure 4: The observed order of accuracy of the code using Eq. 26 for both cases shown in Tab. 4. In first case, the L norm isparallel to the (a) second-order slope and in the second case it is parallel to the (c) fourth-order slope. The resulting curves forthe L norm evolution in time are shown for the (b) first and (d) second case. Both L norm evolutions start from values closeto the machine precision ( − ), since we start the convergence analysis simulations from the exact solution. In the first casethe formal order of accuracy is achieved, p ≈ . , while in the second case p ∗ ≈ . due to discretization error canceling. A truncation error canceling could be seen when the choice for the manufactured solution wavenumberwas q = √ q . The latter modifies the expected convergence rate up to fourth-order, which is consistentwith the analytical development of the manufactured solution in Eq. 23. To verify the general case, weconstruct a manufactured solution with q (cid:54) = √ q . Following this approach, a second-order convergence rateis recovered as theoretically expected.The numerical scheme has been proven to be consistent, unconditionally stable and truly second-orderaccurate in space, in the general case. The chosen grid resolution for the numerical experiments, presentedin the next section, was g r = 16 , which has a good trade-off between resolution and computational cost inorder to represent periodic solutions. 14 .3. Lyapunov functional decay The decay of the “free energy" (Lyapunov) functional is monitored so that we can confirm it is alwaysmonotonically decreasing until a minimum value is reached in the steady state. In order to take into accountthe assumed boundary conditions adopted in this work, we expand the bulk integrals to show that boundaryintegrals vanish, leading to a new expression for the free energy functional, which is going to be discretized.Using integration by parts and the 2D Gauss theorem (first Green identity), we have: (cid:90) Ω d x (cid:0) ψ ∇ ψ (cid:1) = (cid:73) ∂ Ω ψ ∂ψ∂n dl − (cid:90) Ω d x ( ∇ ψ ) , (27)where the first integral from RHS vanishes for the assumed GDBC, above mentioned. Now Eq. 1 can berewritten as: F [ ψ ] = (cid:90) Ω d x (cid:26) − (cid:15) ( x ) ψ + α (cid:2) q ψ − q ( ∇ ψ ) + ( ∇ ψ ) (cid:3) − β ψ + γ ψ (cid:27) . (28)The Lyapunov functional associated to the SH equation is implemented through the discrete formula derivedby Christov & Pontes (2001)[29] for the cubic version, and now extended for the quintic one. The formulapresents a O (cid:0) ∆ t + ∆ x + ∆ y (cid:1) approximation of the functional given by Eq. 2. As pointed by those authors,the monotonic decay of the finite differences version is enforced, provided that the internal iterations converge.The decay is monitored by: F n +1 − F n ∆ t = − n x (cid:88) i =1 n y (cid:88) j =1 (cid:32) ψ n +1 i,j − ψ ni,j ∆ t (cid:33) (29) F n = n x (cid:88) i =1 n y (cid:88) j =1 (cid:20) − (cid:15) (cid:0) ψ ni,j (cid:1) − β (cid:0) ψ ni,j (cid:1) + γ (cid:0) ψ ni,j (cid:1) + αq (cid:0) ψ ni,j (cid:1) (cid:21) − αq n x (cid:88) i =1 n y (cid:88) j =1 (cid:20) ψ ni +1 ,j − ψ ni,j ∆ x (cid:21) + (cid:20) ψ ni,j − ψ ni − ,j ∆ x (cid:21) + (cid:20) ψ ni,j +1 − ψ ni,j ∆ y (cid:21) + (cid:20) ψ ni,j − ψ ni,j − ∆ y (cid:21) + α (cid:88) i =1 h n x n y (cid:88) j =1 (cid:20) ψ ni +1 ,j − ψ ni,j + ψ ni − ,j ∆ x + ψ ni,j +1 − ψ ni,j + ψ ni,j − ∆ y . (cid:21) In order to check if the Lyapunov functional is correctly implemented, we ran two simulations with thecubic SH (with ramped (cid:15) ) and four with the quintic one (both uniformly and ramped (cid:15) ), and monitored theevolution of the functional, of ˙ L and of the patterns observed through the order parameter ψ ( x , t ) . Thepatterns evolution and the associated ˙ L and Lyapunov potentials are presented in Figs. 5, 7, and 8. Thesimulations were run with ramps of the control parameter (cid:15) ( x ) along the x -direction for both GDBC andPBC. The two simulations with the SH3 model started from the same pseudo-random initial condition. Thefour simulations with the SH35 model started from a squared localized patch in the center of the cell. Thispre-existing structure of stripes parallel to the y -direction were constructed by: ψ ( x ,
0) = A cos( q x ) , (30)where A = (cid:110) ( β/ (2 γ ))[1 + (cid:112) (cid:15)γ/β ] (cid:111) / , which is compatible with the expected field ψ amplitude fora stable spatially homogeneous states, for the SH35 [12]. The results of these six simulations confirm thecorrect implementation of the Lyapunov functional for both equations, with a monotonic decay in all cases.15 = 0 t = 10 t = 100 t = 150 t = 200 t = 11885 (a) GBDC t = 0 t = 10 t = 100 t = 150 t = 200 t = 10219 (b) PBC + ˙ L [-] GDBCPBC L y a punov f un c ti on a l ( F [ ψ ]) GDBCPBC (c)
Figure 5: The results of two simulations run to verify the correctness of the implementation of the Lyapunov functional for theSH3 model. First and second rows: SH3 pattern evolution for a ramped system forced with . ≤ (cid:15) ≤ . GDBC (a) and PBC(b), respectively, until the indicated steady state.. Both simulations started with the same pseudo-random initial condition. ˙ L and the Lyapunov potential evolutions are showed below. = 0 t = 10 t = 100 t = 150 t = 200 t = 331 (a) GDBC t = 0 t = 10 t = 100 t = 150 t = 200 t = 328 (b) PBC + ˙ L [-] GDBCPBC L y a punov f un c ti on a l ( F [ ψ ]) GDBCPBC (c)
Figure 6: The results of two simulations run to verify the correctness of the implementation of the Lyapunov functional forthe SH35 model. First and second rows: SH35 pattern evolution for a uniformly forced system with (cid:15) = − . , GDBC (a) andPBC (b), respectively, until the indicated steady state.. Both simulations started from the same pre-existing structure of rollsperpendicular to the gradient of the control parameter. ˙ L and the Lyapunov potential evolutions are showed below. = 0 t = 50 t = 300 t = 1000 t = 5000 t = 15921 (a) GDBC t = 0 t = 50 t = 300 t = 1000 t = 5000 t = 18144 (b) PBC + ˙ L [-] GDBCPBC L y a punov f un c ti on a l ( F [ ψ ]) GDBCPBC (c)
Figure 7: The results of two simulations run to verify the correctness of the implementation of the Lyapunov functional for theSH35 model. First and second rows: SH35 pattern evolution for a ramped system forced with − . ≤ (cid:15) ≤ − . , GDBC (a) andPBC (b), respectively, until the indicated steady state.. Both simulations started from the same pre-existing structure of rollsperpendicular to the gradient of the control parameter. ˙ L and the Lyapunov potential evolutions are showed below. = 0 t = 10 t = 1000 t = 1500 t = 2500 t = 11378 (a) GDBC t = 0 t = 10 t = 100 t = 170 t = 260 t = 22199 (b) PBC + ˙ L [-] GDBCPBC L y a punov f un c ti on a l ( F [ ψ ]) GDBCPBC (c)
Figure 8: The results of two simulations run to verify the correctness of the implementation of the Lyapunov functional for theSH35 model. First and second rows: SH35 pattern evolution for a ramped system forced with − . ≤ (cid:15) ≤ − . , GDBC (a) andPBC (b), respectively, until the indicated steady state.. Both simulations started from the same pre-existing structure of rollsperpendicular to the gradient of the control parameter. ˙ L and the Lyapunov potential evolutions are showed below. . Numerical experiments The numerical scheme described in Sec. 2 was used to solve the SH3 and the SH35 equations in squaredomains with rigid and periodic boundary conditions, and with nonuniform forcings. The parameter valuesadopted are given in Tabs. 1, 5 and 6. All simulations were done in a satisfactory grid resolution of thespatial grid of points per critical wavelength ( g r = 16 ). In the case of the SH3 equation we present the result of ten simulations (denoted as configurations to ),five of them with GDBC and five with PBC. The resulting steady state for each simulation are summarizedin Fig. 9. All simulations started with the same pseudo-random initial condition. For each type of boundarycondition we ran a configuration with uniform forcing to reproduce existing results, two configurations withramps of the control parameter (cid:15) , as detailed in Tab. 5, and two ones with different gaussian distributions ofthe control parameter. In ramped systems, the parameter is assumed as uniform along the y -direction. Thegaussian distribution of (cid:15) introduces symmetric gradients in the radial direction.The steady state solutions obtained are shown in Fig. 9. The first line in this figure shows the distributionof the control parameter (cid:15) along the x -direction of the domain. Second and fourth lines present the steadystate solutions with rigid and with periodic boundary conditions, respectively. Third an fifth lines presentthe same results shown in lines two and four, but this time with frames constructed with enhanced contrastto make visible the subcritical structure developed in regions where the control parameter (cid:15) is negative. Eachconfiguration run is identified with an assigned number.Configuration 1 shows a pattern emerging in a uniformly forced system with GDBC, in qualitative agree-ment with previous results [22, 25, 26, 28]. The case was run in the verification framework our numericalcode. A pattern with unavoidable defects is developed, in order to match the strong requirement of rollsapproaching perpendicularly the side walls. A structure with lower density of defects appears in configuration6, where the same forcing of Frames 1 is assumed, but now, with periodic boundary conditions. The patterndevelops a zig-zag instability.Configuration 2 reproduces the result of a pattern developed in presence of a ramped control parameter (cid:15) , with a subcritical region [27, 29, 28]. A structure with smaller amount of defects emerges, thanks to thefact that the weak structure of rolls parallel to the wall appears in the subcritical region. A similar situationoccurs in configuration 3, where the system is forced with a ramp of the control parameter (cid:15) , however takingnon negative values. In this case an weak structure of rolls perpendicular to the lower ( y = 0) sidewall isvisible.Configurations 7 and 8 present the results for systems with same forcing of configurations 2 and 3, butwith periodic boundary conditions. A tendency to develop a structure of rolls parallel to the gradient of thecontrol parameter (cid:15) clearly appears, as well as the existence of a Benjamin-Feir instability close to the leftwall, induced by the the supercritical region close to the right wall.20 able 5: Parameters assumed for the numerical experiments with the SH3 equation. Parameter Formulæ Value Description q - 1.0 Critical wavenumber λ π/q π Critical wavelength w x , w y - 10 Wavelengths per domain length g r - 16 Grid resolution (nodes per wavelength) n x , n y w x × g r , w y × g r
160 Nodes per mesh side ( n ) N n x × n y × Total number of mesh nodes L x , L y w x λ , w y λ ≈ . Square domain length ( L ) ∆ x , ∆ y L/ ( n − ≈ . Space step (GDBC) ∆ x , ∆ y L/n ≈ . Space step (PBC) ∆ t - 0.5 Time step for SH3 A - 0.2 Gaussian maximum value (peak) R N - For configurations 4 and 9 R . N - For configurations 5 and 10 x , y L x / , L y / - Gaussian centerConfigurations 4, 5, 9 and 10 present the result of simulations performed with a Gaussian radial distri-bution of the control parameter (cid:15) , in the form: (cid:15) ( x ) = Ae − R (( x − x ) +( y − y ) ) , (31)where the parameters are specified in Tab. 5. The parameter R is assumed to be R in configurations 4 and9; and R in configurations 5 and 10. Worth noticing that these gaussian distributions introduce a gradientof (cid:15) in all directions with the maximum of (cid:15) at the center of the domain.Configurations 4 and 5 were run with rigid boundary conditions and configurations 9 and 10, with periodicconditions. A sharp Gaussian, rapidly decreasing from the maximum value across a short distance, leadsto the onset of a target pattern with both boundary conditions. The target pattern collapses in presenceof a wider Gaussian distribution of the control parameter (cid:15) , leading to a structure of rolls. The structurewavevector takes the direction of one of the domain diagonals, again, since the most difficult direction formodulations of the structure is the one associated to the wavevector, whereas the easiest one is the directionperpendicular to the wavevector.The time evolution of configurations 3 (GDBC) and 8 (PBC), from the pseudo-random initial conditionto the steady state in shown in Fig. 10. The evolution of configurations 5 (GDBC)and 10 (PBC) are shownin Fig. 11. All curves show the expected result, with ˙ L rapidly decreasing as a pattern emerges and the21mplitude of the structure essentially saturates. This stage is followed by a much slower phase dynamics,ending with an exponential decay towards the steady state. In general, ˙ L peaks are associated to the removalof defects. We observe that boundary conditions have little effect on the resulting curves. x − . . . (cid:15) x − . . . (cid:15) x − . . . (cid:15) x , y − . . . (cid:15) x , y − . . . (cid:15) t = 105262 t = 27252 t = 39972 t = 5780 t = 319880 t = 5560 t = 7270 t = 5377 t = 8325 t = 84100 Figure 9: Zoology of ten patterns obtained by numerical integration of the SH3 equation with GDBC (rows 1 and 2) andPBC (rows 3 and 4). All patterns shown correspond to the steady attained when ˙ L ≤ . × − (Eq. 21). All simulationsstarted from the same pseudo-random initial condition. The applied forcing is shown in the first row. First column: uniformforcing; Second and third columns: ramps of the control parameter (cid:15) ; Fifth and sixth columns: Two Gaussian forcings, accordingto Eq.31. The first and third rows correspond to GDBC and PBC, respectively. Second and fourth rows contains the samepatterns shown in rows one and three, respectively, with an enhanced contrast. Patterns configurations 2 and 3 tend to approachperpendicularly to the supercritical boundaries and parallel to subcritical ones. The pattern configuration 4 presents a targetstructure compatible with the narrow Gaussian profile of the forcing. The target collapses in the pattern configuration 5, andevolves to a structure of rolls parallel to one of diagonals, and not to the side walls. This effect occurs as a consequenceof the modulation, which if harder along the direction of the pattern wavevector. Similar phenomena occurs with patternsconfigurations 9 and 10. = 0 t = 50 t = 100 t = 250 t = 350 t = 350 t = 380 t = 400 t = 410 t = 500 t = 750 t = 37500 (a) GDBC t = 0 t = 50 t = 100 t = 300 t = 500 t = 750 t = 1000 t = 1250 t = 1500 t = 1750 t = 2000 t = 5377 (b) PBC ˙ L [-] GDBCPBC (c)Figure 10: First and second lines: pattern evolution for the configuration 3 of Fig. 9, with GDBC (a), until the steady state.Third and fourth lines: pattern evolution for the configuration 8 of Fig. 9, with PBC (b), until the steady state. In both cases,the system is forced with a ramp of the control parameter (cid:15) along the x direction, with − . ≤ (cid:15) ≤ . . The associated ˙ L × t curves are also shown (c). = 0 t = 200 t = 27700 t = 27850 t = 28000 t = 28010 t = 28040 t = 28050 t = 28060 t = 28080 t = 28100 t = 319880 (a) GDBC t = 0 t = 10 t = 20 t = 30 t = 40 t = 50 t = 100 t = 200 t = 300 t = 400 t = 2000 t = 8325 (b) PBC ˙ L [-] GDBCPBC (c)Figure 11: First and second lines: pattern evolution for the configuration 5 of Fig. 9, with GDBC (a), until the steady state.Third and fourth lines: pattern evolution for the configuration 10 of Fig. 9, with PBC (b), until the steady state. In both cases,the system is forced with a gaussian distribution of the control parameter (cid:15) . The associated ˙ L × t curves are also shown (c). .2. Numerical experiments with the SH35 equation The quintic version of the Swift-Hohenberg equation with destabilizing cubic term allows for the existenceof stable localized solutions for an interval of subcritical values of (cid:15) , i.e. coexistence of the solution ψ = 0 and a modulated one, shown by [12]. The parameter values adopted are given in Tab.6. In this sense,we run six experiments with the SH35 equation: two with uniform forcing ( (cid:15) = − . , GDBC and PBC,respectively, two with a ramp − . ≤ (cid:15) ≤ − . , GDBC and PBC, respectively, and two last ones with a ramp − . ≤ (cid:15) ≤ − . , GDBC and PBC, respectively. The experiments with uniform forcing were conducted ascode verification, to reproduce existing results. The time evolution of these two experiments are shown inFig. 6. Both simulations started from the same pseudo-random initial condition.The four experiments with ramped systems started from the same pre-existing structure of rolls perpen-dicular to the gradient of the control parameter. These experiments were used to verify the implementationof the Lyapunov functional with the quintic equation SH35. The evolution of these simulations are shown inFig. 7, and 8. In the case of a ramp given by − . ≤ (cid:15) ≤ − . the pre-existing structure evolves towards alocalized structure of rolls with the same orientation of the initial condition, with both boundary conditions.When the applied forcing takes the form − . ≤ (cid:15) ≤ − . , lowering the energy associated with the non-trivialsolution both with GDBC and PBC, a structure of rolls with the same orientation of the initial conditionalso prevails, but the final structure is localized with GDBC and occupies the entire domain, with PBC. Table 6: Parameters assumed for the numerical experiments with the SH35 equation.
Parameter Formulæ Value Description q - 1.0 Critical wavenumber λ π/q π Critical wavelength w x , w y - 8 Wavelengths per domain length g r - 16 Grid resolution (nodes per wavelength) n x , n y w x × g r , w y × g r
128 Nodes per mesh side ( n ) N n x × n y × Total number of mesh nodes L x , L y w x λ , w y λ ≈ . Square domain length ( L ) ∆ x , ∆ y L/ ( n − ≈ . Space step (GDBC) ∆ x , ∆ y L/n ≈ . Space step (PBC) ∆ t - 0.1 Time step for SH3525 − . − . − . (cid:15) x − . − . (cid:15) x − . − . (cid:15) t = 331 t = 15921 t = 11378 t = 328 t = 18144 t = 22199 Figure 12: Zoology of six patterns obtained by numerical integration of the SH35 equation with GDBC (rows 1 and 2) and PBC(rows 3 and 4). All patterns shown correspond to the steady attained when ˙ L ≤ . × − (Eq. 21) All simulations startedfrom a pre-existing structure of stripes oriented along the y direction. The applied forcing is shown in the first row. First column:uniform forcing with (cid:15) = − . . Second column: forcing with ramp of the control parameter in the form − . ≤ (cid:15) ≤ − . . Thirdcolumn: ramp in the form − . ≤ (cid:15) ≤ − . . The first and third rows correspond to GDBC and PBC, respectively. Secondand fourth rows contain the same patterns shown in rows one and three, respectively, with an enhanced contrast. In all cases,the pre-existing structure of stripes oriented along the y prevails at the steady state. However, configurations 11, 12, 14 and 15evolved to localized structures, whereas, in the case of configurations 13 and 16 the resulting pattern occupy the entire domain. In the present endeavor, we extended a numerical scheme proposed by Christov & Pontes [28] to investigatepattern formation modelled by the cubic Swift-Hohenberg equation in two dimensions. The original schemepresents second order representation of all derivatives, strict implementation of the associated Lyapunovfunctional, rigid boundary conditions (GDBC), and a semi-implicit time discretization of the terms, aiming26t constructing negative definite operators that assure the unconditional stability of the scheme.The present work includes the quintic version of the Swift-Hohenberg equation and PBC for both thecubic and the quintic versions of the model. The scheme retains all characteritics of the original one, namelystrict representation of the Lyapunov functional, unconditional stability, and second order representation ofall derivatives. In addition, we also included a convergence analysis of the scheme, verification tests, andan initial evaluation of the effect of nonuniform forcings in the resulting pattern. The nonuniform forcingsadopted in the present work consist of spatial ramps and of gaussian distributions of the control parameter (cid:15) . A detailed study of the effects of nonuniform forcings will be addressed in a forthcoming paper.Verification tests consisted in qualitatively reproducing existing results of simulations with the SH3 andSH35 equations with uniform forcings and both GDBC and PBC. Good agreement was achieved. As anadditional verification procedure, we also found good agreement between the results with the SH3 forcedwith a spatial ramp of (cid:15) , and those presented by Pontes et al. (2008)[28]. In that work, the authors adoptedGDBC and a first order in time scheme. The results of the SH3 equation with GDBC and ramps of (cid:15) arepresented in Figs. 2, 5a and in simulations configurations 1, 2 and 3 of Fig. 12.Rigid boundary conditions (GDBC) make patterns of supercritical stripes developed close to differentboundaries to compete, due to the tendency of stripes to aproach boundaries perpendicularly to them.Stripes cannot be perpendicular to all boundaries simultaneously. In addition, boundary effects competewith bulk ones. The overall result is a structures with high density of deffects. A different situation occurswhen PBC are imposed. In this case, a periodicity is imposed to the emerging pattern, and boundary effectsare excluded. The result is, in general, a pattern with fewer defects than ones developed with the sameforcing and GDBC. This is also the case when the system is forced with nonuniform distributions of thecontrol parameter.An important configuration in nonuniformly forced systems consists in adopting GDBC, and a forcingwhere part of the domain is mantained at subcritical conditions. In this case a less known result points to thefact that the subcritical pattern of stripes, induced by the supercritical region, approach sidewalls parallelto them. These patterns are automatically perpendicular to adjacent walls. As a result, nonuniformlyforced systems with a subcritical region tend to develop paterns with a lower density of defects than systemsuniformly forced [28]. The presented simulations reproduce this effect.A particular feature observed in nonuniformly forced systems, using the SH3 equation, PBC and spatialramps of the control parameter, consists in the development of several examples of the Benjamin-Feir insta-bility, which manifest itself in the form of “alternating” rolls, in a dislocation sense, promoted by a regionwhere minimum and maximum values of (cid:15) are close.The quintic SH35 equation with a destabilizing cubic term allows for the emergence of localized solutionsin an interval of subcritical values of (cid:15) . Six simulations were run with the quintic equation, assuming threedifferent forcings and two boundary conditions for each one. All simulations started from the same pseudo-random initial condition. The resulting steady states are presented in Fig. 12. The first group consisted27f two simulations with uniform forcing (cid:15) = − . . These simulations are denoted by configurations 11 and14 in Fig. 12 The second group was forced with a ramp along the x direction, given by − . ≤ (cid:15) ≤ . (configurations 12 and 15). The third group was run with . ≤ (cid:15) ≤ . (configurations 13 and 16). Inall configurations, the resulting patterns consisted of stripes perpendicular to the gradient of the controlparameter, with localized patterns in the first and second groups, and the stripes occupying the entiredomain in the third group.A mesh of g r = 16 points per critical wavelength was adopted in all simulations, which represents a goodtrade-off between spatial resolution and computational cost.
5. Conclusion
In this article, we extended a numerical scheme proposed by Christov & Pontes [28] to investigate patternformation modeled by the cubic Swift-Hohenberg equation SH3 in two dimensions. The original schemepresents second order representation of all derivatives, strict implementation of the associated Lyapunovfunctional, rigid boundary conditions (GDBC), and a semi-implicit assignment of the terms.The present work includes the quintic version of the Swift-Hohenberg equation and periodic boundaryconditions (PBC) for both the cubic and the quintic versions of the model. The scheme retains all charac-teristics of the original one, namely strict representation of the Lyapunov functional, unconditional stability,and second order representation of all derivatives. In addition, we also included a convergence analysis, newverification tests, and an initial evaluation of the effect of nonuniform forcings in the form of spatial rampsand of gaussian distributions of the control parameter (cid:15) . Among the verification tests, the convergence anal-ysis confirmed the truly second-order accuracy of the scheme both in space and time and the existence oflocalized structures developed in the framework of the SH35 equation even with nonuniform forcings. Thenumerical experiments conducted in this work suggest the existence of effects and the onset of patterns notaddressed in the literature. Both questions will be the object of a forthcoming paper.The conducted tests confirmed the robustness of the developed tool for pursuing the investigation ofthe pattern formation through the parabolic Swift-Hohenberg equation presenting various nonlinear terms,including nonuniform forcings.
6. Acknowledgments
The authors thank FAPERJ (Research Support Foundation of the State of Rio de Janeiro) and CNPq(National Council for Scientific and Technological Development) for the financial support. Daniel Coelhoacknowledges a fellowship from the Coordination for the Improvement of Higher Education Personnel-CAPES(Brazil). A FAPERJ Senior Researcher Fellowship is acknowledged by J. Pontes. The authors dedicate aspecial posthumous thanks to C. I. Christov, who performed great development regarding the numericalmethods addressed in this paper. 28 eferences [1] L. Yang, A. M. Zhabotinsky, I. R. Epstein, Stable squares and other oscillatory turing patterns in areaction-diffusion model, Physical review letters 92 (2004) 198303.[2] J. V. Moloney, A. C. Newell, Nonlinear Optics, Westview Press, 2008.[3] D. Walgraef, Spatio-Temporal Pattern Formation: With Examples from Physics, Chemistry, and Mate-rials Science, Partially Ordered Systems, 1 ed., Springer-Verlag New York, 1997.[4] N. Ghoniem, D. Walgraef, Instabilities and Self-organization in Materials, Oxford Univ. Press, 2008.[5] N. Provatas, K. Elder, Phase-Field Methods in Materials Science and Engineering, Wiley, 2010.[6] M. Plapp, A. Karma, Multiscale finite-difference-diffusion-monte-carlo method for simulating dendriticsolidification, Journal of Computational Physics 165 (2000) 592–619.[7] N. Provatas, M. Greenwood, B. Athreya, N. Goldenfeld, J. Dantzig, Multiscale modeling of solidification:phase-field methods to adaptive mesh refinement, International Journal of Modern Physics B 19 (2005)4525–4565.[8] S. Ghosh, M. Plapp, Influence of interphase boundary anisotropy on bulk eutectic solidification mi-crostructures, Acta Materialia 140 (2017) 140–148.[9] S. Ghosh, A. Karma, M. Plapp, S. Akamatsu, S. Bottin-Rousseau, G. Faivre, Influence of morphologicalinstability on grain boundary trajectory during directional solidification, Acta Materialia 175 (2019)214–221.[10] J. Swift, P. C. Hohenberg, Hydrodynamic fluctuations at the convective instability, Physical Review A15 (1977).[11] J. Pontes, D. Walgraef, E. Aifantis, On dislocation patterning: Multiple slip effects in the rate equationapproach, International journal of plasticity 22 (2006) 1486–1505.[12] H. Sakaguchi, H. R. Brand, Stable localized solutions of arbitrary length for the quintic Swift-Hohenbergequation, Physica D 97 (1996) 274–285.[13] I. Walton, On the onset of Rayleigh-Bénard convection in a fluid layer of slowly increasing depth, Studiesin Applied Mathematics 67 (1982) 199–216.[14] I. Walton, The onset of cellular convection in a shallow two-dimensional container of fluid heatednon-uniformly from below, Journal of Fluid Mechanics 131 (1983) 455–470.[15] K. Elder, M. Katakowski, M. Haataja, M. Grant, Modeling elasticity in crystal growth, Physical reviewletters 88 (2002) 245701. 2916] K. Elder, M. Grant, Modeling elastic and plastic deformations in nonequilibrium processing using phasefield crystals, Physical Review E 70 (2004) 051605.[17] J. Langer, H. Müller-Krumbhaar, Theory of dendritic growth—i. elements of a stability analysis, ActaMetallurgica 26 (1978) 1681–1687.[18] J. S. Langer, Instabilities and pattern formation in crystal growth, Reviews of modern physics 52 (1980)1.[19] M. Elsey, B. Wirth, A simple and efficient scheme for phase field crystal simulation, ESAIM: Mathe-matical Modelling and Numerical Analysis 47 (2013) 1413–1432.[20] H. G. Lee, J. Kim, A simple and efficient finite difference method for the phase-field crystal equationon curved surfaces, Computer Methods in Applied Mechanics and Engineering 307 (2016) 32–43.[21] Y. Li, J. Kim, An efficient and stable compact fourth-order finite difference scheme for the phase fieldcrystal equation, Computer Methods in Applied Mechanics and Engineering 319 (2017) 194–216.[22] H. S. Greenside, J. W. M. Coughran, Nonlinear pattern formation near the onset of Rayleigh-Bénardconvection, Physical Review A 30 (1984) 398–428.[23] M. Cross, Ingredients of a theory of convective textures close to onset, Physical Review A 25 (1982)1065.[24] M. Cross, Boundary conditions on the envelope function of convective rolls close to onset, The Physicsof Fluids 25 (1982) 936–941.[25] P. Manneville, Non-equilibrium Structures and Weak Turbulence, Academic Press: San Diego, CA, USA,1990.[26] M. C. Cross, Pattern formation outside equilibrium, Rev. Modern Physics 65 (1983) 851–1112.[27] C. Christov, J. Pontes, D. Walgraef, M. Velarde, Implicit time splitting for fourth-order parabolicequations, Computer Methods in Applied Mechanics and Engineering 148 (1997).[28] J. Pontes, D. Walgraef, C. I. Christov, Pattern formation in spatially ramped Rayleigh-Bénard systems,Journal of Computational Interdisciplinary Sciences 1 (2008) 11–32.[29] C. Christov, J. Pontes, Numerical scheme for Swift-Hohenberg equation with strict implementation ofLyapunov functional, Mathematical and Computer Modelling 35 (2002).[30] H. G. Lee, An energy stable method for the swift–hohenberg equation with quadratic–cubic nonlinearity,Computer Methods in Applied Mechanics and Engineering 343 (2019) 40–51.3031] E. Vitral, P. H. Leo, J. Viñals, Role of gaussian curvature on local equilibrium and dynamics of smectic-isotropic interfaces, Physical Review E 100 (2019) 032805.[32] J. Wang, S. Zhai, A fast and efficient numerical algorithm for the nonlocal conservative Swift–Hohenbergequation, Mathematical Problems in Engineering (2020).[33] J. Douglas, H. H. Rachford., On the numerical solution of heat conduction problems in two and threespace variables, Trans. Amer. Math. Soc. (1956).[34] N. N. Yanenko, M. Holt, The Method of Fractional Steps, Springer, 1971.[35] E. Vitral, D. Walgraef, J. Pontes, G. Anjos, N. Mangiavacchi, Nano-patterning of surfaces by ionsputtering: Numerical study of the anisotropic damped Kuramoto-Sivashinsky equation, ComputationalMaterials Science 146 (2018).[36] P. J. Roache, Code verification by the method of manufactured solutions, J. Fluids Eng. 124 (2002)4–10.[37] C. J. Roy, Review of code and solution verification procedures for computational simulation, Journal ofComputational Physics 205 (2005) 131–156.[38] V. E. Podneks, I. W. Hamley, Landau-Brazovskii theory for the Ia ¯3 dd