Finite element solution of nonlocal Cahn-Hilliard equations with feedback control time step size adaptivity
Gabriel F. Barros, Adriano M. A. Côrtes, Alvaro L. G. A. Coutinho
FF INITE ELEMENT SOLUTION OF NONLOCAL C AHN -H ILLIARDEQUATIONS WITH FEEDBACK CONTROL TIME STEP SIZEADAPTIVITY
A P
REPRINT
Gabriel F. Barros
COPPEFederal University of Rio de JaneiroRio de Janeiro, RJ 21941-598 Brazil [email protected]
Adriano M. A. Cortes
NUMPEX-COMPFederal University of Rio de Janeiro, BrazilDuque de Caxias, RJ 25240-005 Brazil [email protected]
Alvaro L. G. A. Coutinho
COPPEFederal University of Rio de JaneiroRio de Janeiro, RJ 21941-598 Brazil [email protected]
October 1, 2020 A BSTRACT
In this study, we evaluate the performance of feedback control-based time step adaptivity schemes forthe nonlocal Cahn-Hilliard equation derived from the Ohta-Kawasaki free energy functional. Thetemporal adaptivity scheme is recast under the linear feedback control theory equipped with an errorestimation that extrapolates the solution obtained from an energy-stable, fully implicit time marchingscheme. We test three time step controllers with different properties: a simple Integral controller, acomplete Proportional-Integral-Derivative controller, and the PC11 predictive controller. We assessthe performance of the adaptive schemes for the nonlocal Cahn-Hilliard equation in terms of thenumber of time steps required for the complete simulation and the computational effort measured bythe required number of nonlinear and linear solver iterations. We also present numerical evidence ofmass conservation and free energy decay for simulations with the three different time step controllers.The PC11 predictive controller is the best in all three-dimensional test cases. K eywords Nonlocal Cahn-Hilliard equation · Time step size adaptivity · Ohta-Kawasaki Functional · FeedbackControl Theory
The Cahn-Hilliard equation (or simply the CH equation) was derived in 1958 to model the phase separation in binaryalloys [1, 2]. Since then the CH equation appears in several different physical contexts [3] such as diblock copolymers[4, 5], image inpainting [6], binary fluid flow [7, 8, 9, 10, 11], fracture propagation [12, 13], tumour growth [14, 15]and topology optimization [16], to mention a few. The CH equation is ∂φ∂t = ∇ · (cid:20) M ( φ ) ∇ (cid:18) ∂ Ψ ∂φ − (cid:15) ∆ φ (cid:19)(cid:21) , (1)where φ ( x , t ) is a field evolving in time and space, M ( φ ) > is the mobility, related to the diffusion process, Ψ( φ ) isthe bulk free energy, and (cid:15) is a parameter related to the interfacial energy. From a strict mathematical point of view,considering no-flux or periodic boundary conditions, the equation can also be seen as the gradient flow [17] of the a r X i v : . [ c s . C E ] S e p PREPRINT - O
CTOBER
1, 2020Ginzburg-Landau functional described by F ( φ ) = (cid:90) Ω (cid:18) Ψ( φ ) + (cid:15) |∇ φ | (cid:19) d Ω . (2)The derivation of (1) considers only short-ranged microforces. Therefore, the use of the CH equation becomes restrictedto model physical phenomena where only local interactions of particles are taken into account. The development of anonlocal Cahn-Hilliard equation (NCH for short) [18, 19, 20] fills this gap through the derivation of a phase-separationmodel that also considers long-range interactions. Among the various NCH applications, we highlight the modeling ofthe diblock copolymer self-assembly. Diblock copolymers are a specific class of copolymers, in which two chemicallydistinct monomer units are grouped in discrete blocks along the polymer chain [21]. The large variety of morphologiesobtained in the resulting self-assembled copolymer by manipulating different molecular parameters, added to thecrescent research interest in nanotechnology, lead to the development of block copolymers related discoveries inadvanced materials, drug delivery, patterning, porous materials, and many others over the last decades [22]. Sincewe are interested in the patterns formed by the solution of the NCH equation, one of the key questions in the diblockcopolymer context that remains open is if it is possible to find which pattern minimizes the energy over all possiblepatterns for a given set of parameters. The Ohta-Kawasaki (O–K) free energy functional models pattern morphologiesvia energy minimization involving the competition of both short and long-range microforces. In [23, 24, 25], wefind analytical solutions in one dimension for a global minimizer for the O–K free energy functional. However, theyare of limited use. Numerical simulations of the NCH equation for different parameter sets became fundamental inexploring the associated phase diagrams for diblock copolymer melts [5, 26, 27, 28, 29]. However, the use of thesenumerical simulations is nontrivial and can present difficulties related to the stiffness of the equations and demand largecomputational power [30]. In this sense, the development of computational techniques that improve the performanceand accuracy of the computational modeling of the NCH equation is an active topic of research.This paper evaluates time step adaptive schemes for the NCH equation in the linear feedback control theory context withproper error estimation and time step controller. Several works presented adaptive time-stepping using error estimationand time step controllers in the CH equation [31, 32, 33, 34, 35], the use of more sophisticated controllers and errorestimation techniques are not widely addressed in the nonlocal case. In this study, we use three different time-stepcontrollers with different properties and behavior, together with a proper error estimation method that prevents thecalculation of the same time step multiple times. We evaluate the controllers in terms of performance in differentexamples in terms of required time steps for the completion of the simulations and the total number of nonlinearand linear iterations required by the solvers. We also assess the employment of an implicit time-marching scheme,mathematically and numerically proven to be energy-stable in standard phase-field functionals. Our results revealnumerical evidence of mass conservation and free energy decay for the nonlocal case as anticipated theoretically in[5, 36, 37]. The paper is structured as follows. Section 2 describes the NCH formulation. Numerical and computationalimplementation details are shown in Section 3. Section 4 introduces the temporal adaptivity schemes, and in Section 5,we show several numerical examples for solving the CH/NCH equations with adaptive time step control. The paperends with the conclusions we drew from the comparison of the three controllers in different contexts. We consider the following nonlocal extension of the standard CH equation, that is, the NCH equation, ∂φ∂t = ∇ · (cid:20) M ( φ ) ∇ (cid:18) ∂ Ψ ∂φ − (cid:15) ∆ φ (cid:19)(cid:21) − σ ( φ − ¯ φ ) , (3)where σ represents the nonlocal parameter, responsible for modeling the magnitude of the long-range microforcesbetween the phases. The parameter ¯ φ is the mean value of φ in the domain Ω ∈ R n sd with boundary ∂ Ω and n sd = 2 , ,such that, ¯ φ = 1 | Ω | (cid:90) Ω φd Ω . (4)One way of deriving Eq. (3) is under the H − (Ω) gradient flow from the O–K free energy functional [38]. The O–Kfunctional can be written [27, 5] as, F ( φ ) = (cid:90) Ω (cid:18) Ψ( φ ) + (cid:15) |∇ φ | (cid:19) d Ω + (cid:90) Ω σ |∇ v | d Ω , (5)where v is related to φ via the boundary value problem − ∆ v = ( φ − ¯ φ ) .2 PREPRINT - O
CTOBER
1, 2020The O–K free energy functional is derived from the mean field theory in the context of diblock copolymers [4]. In thecase where σ = 0 , the O–K functional (5) becomes the Ginzburg-Landau free energy functional (2). Consequently,in this case, the NCH equation (3) becomes the CH equation (1). Both CH [34] and NCH [39] equations minimizethe interfaces between phases, solving the isoperimetric problem. The difference, however, lies in the fact that theGinzburg-Landau free energy functional is minimized through the separation of phases due to short-range microforceswhile the O–K free energy functional models pattern morphologies via energy minimization involving the competitionof both short and long-range microforces, where the latter is modeled by the magnitude of the nonlocal parameter σ . Thecompetition between local and nonlocal microforces in the O–K free energy functional leads to many different patternformations in the equilibrium configuration of a copolymer melt, such as lamellae, spheres, gyroids, and cylinders [21].It is possible to map a given copolymer structure to a set of NCH parameters, such as (cid:15) , ¯ φ , and σ by a phase diagram[5, 27, 28].Equation (3) is solved on a bounded domain Ω with Lipschitz-continuous boundaries ∂ Ω , and on the time interval [0 , T ] with prescribed initial conditions φ ( x ,
0) = φ . Regarding the boundary conditions the usual ones are theno-flux boundary conditions, that is, ∇ φ · n = 0 and M ( φ ) ∇ ( ∂ Ψ ∂φ − (cid:15) ∇ φ ) · n = 0 , and the periodic boundaryconditions. The use of these boundary conditions implies on mass conservation and free energy decay for both CH[40] and NCH [36, 37] equations. In the past decades, several studies investigated numerical strategies to solve theCH equation. Considering spatial discretization, there are models based on finite differences [41], finite volumes[31], finite elements [42, 43, 40] and spectral methods [44], all with their advantages and drawbacks. In the finiteelement method, the presence of fourth-order spatial derivatives in the CH equation (in the strong form) requires theuse of C -continuous elements in the primal variational formulation of the equation. Stogner et al. [35] presented afinite element formulation using C -continuous elements in two dimensions for rectangular grids, while other studiescircumvented this situation with different techniques such as NURBS-based isogeometric analysis, [32], variablesplitting technique (mixed formulation) [43] and discontinuous Galerkin methods [45]. The CH equation temporaldiscretization is also nontrivial. Besides being a stiff and nonlinear equation, which makes it practically unsolvableby explicit methods [41], the time integration method must generally obey an energy decay property since, in mostcases, e.g., considering no-flux or periodic boundary conditions, the free energy functional is a Lyapunov functional.Many studies developed different time integration methods for the CH equation preserving this property [33, 32]. Interms of numerical strategies explicitly designed for the NCH equation, the literature is not as rich as the case of the CHequation. We highlight the development of preconditioners [38, 46, 30], computational implementation details, andtheoretical proofs regarding stability, boundedness, and mass conservation using the finite element method [36].Several strategies have been proposed to increase the accuracy and performance of the simulations. For example, thespatial discretization of the CH equation must be fine enough to consider the smooth transition of the interface thatarises between different phases [34] while the bulk domain does not require fine meshes. Therefore, it is commonto track the interface areas to refine the mesh while coarsening the bulk domain [31]. When considering temporaldiscretization, some physical phenomena described by the CH equation require small time step sizes to capture fastdynamics. However, there are stages where the dynamics are slow, and consequently, the use of larger time stepsis allowed. The use of smaller time step sizes in these stages is translated into unnecessary computational costs.Nevertheless, an adaptive time-stepping strategy can help determine whether or not the time step size can be enlargedor reduced. The same problems also appear in the nonlocal case and, although several studies present strategies tocircumvent these difficulties in the CH context, this issue is not widely addressed in the NCH literature. The presentpaper contributes to fill this gap. In this study, the finite element method is employed to discretize in space the NCH equation. We use for temporalintegration, an implicit, second-order, unconditionally energy-stable method originally proposed for the CH equationand other traditional phase-field equations [33]. This method enables the use of larger time steps obtained by timeadaptivity without affecting the numerical stability of the simulations. All the numerical solutions are computed usingthe
FEniCS framework version 2019.1.0 [47, 48], a high-performance finite element library written in Python/C++.
The formulation of the NCH equation contains a biharmonic operator. Thus standard C -continuous finite elements arenot suitable for its primal variational formulation. Nevertheless, a splitting strategy, also called mixed formulation, isemployed to avoid the continuity constraint and enable the use of C -continuous elements to approximate the solutionof the CH [43] and the NCH [36, 30] equations, converting the nonlinear equation into a coupled nonlinear system,with two degrees of freedom per node. 3 PREPRINT - O
CTOBER
1, 2020The split form can be achieved by introducing the chemical potential µ as an unknown field. Given a spatial domain Ω ∈ R n sd with boundaries ∂ Ω and n sd = 2 , and the primal NCH equation on (1), its split version is given by ∂φ∂t = ∇ · ( M ( φ ) ∇ µ ) − σ ( φ − ¯ φ ) , (6) µ = ∂ Ψ ∂φ − (cid:15) ∇ φ. The weak form of the system can be obtained by integrating both equations (6) in their strong form against weightingfunctions q, w ∈ H (Ω) , where H (Ω) is the Sobolev space of the square integrable functions with an integrable firstweak derivative, and applying the divergence theorem. The Galerkin method approximates the unknown fields throughfunctions in a finite dimension space. Considering a partition of the form Ω = (cid:83) e Ω e , and being P k (Ω e ) the space ofpolynomials of degree equal or less than k over Ω e , the function spaces are defined as: S ht = { φ h ( · , t ) , µ h ( · , t ) ∈ H (Ω) | φ h ( · , t ) | Ω e , µ h ( · , t ) | Ω e ∈ P k (Ω e ) , ∀ e } , (7) W h = { w h , q h ∈ H (Ω) | w h | Ω e , q h | Ω e ∈ P k (Ω e ) , ∀ e } . (8)For a standard finite element discretization, the semi-discrete finite element formulation for the NCH nonlinear systemis: Given φ ( x ,
0) = φ h ( x ) , find φ h ( t ) , µ h ( t ) ∈ S ht , ∀ w h , q h ∈ W h , so that: (cid:18) w h , ∂φ h ∂t (cid:19) + ( ∇ w h , M ( φ h ) ∇ µ h ) + ( w h , σφ h ) − ( w h , σ ¯ φ h ) = 0 , (9) ( q h , µ h ) − (cid:18) q h , ∂ Ψ ∂φ (cid:19) − ( ∇ q h , (cid:15) ∇ φ h ) = 0 , where ( · , · ) is the L inner product. After splitting the NCH equation, the chemical potential µ becomes another solvablefield. Periodic boundary conditions as well as no-flux boundary conditions are considered. The NCH equation is a time-dependent equation, so a proper time integration method must be chosen. It is essential touse an energy stable integration method since both the Ginzburg Landau and the Ohta-Kawasaki free energy functionalsare Lyapunov functionals when no-flux or periodic boundary conditions are applied to the domain [49, 5, 37].The choice of a time integration method for the CH/NCH equations is not a trivial task. Explicit methods are oftenprohibitive due to severe restrictions on the time step size, which is around O (∆ x ) , arising from the stiffness ofthe equations. Fully implicit methods, because of the nonlinear nature of the CH/NCH equations, require nonlinearsolvers, increasing memory and computational requirements. Implicit methods allow larger time steps, but if thetime step is too large, the nonlinear discrete systems emanating from the CH equation can present multiple solutions[41, 49]. An intermediate approach is provided by semi-implicit time-stepping algorithms, where some terms areimplicitly treated while others remain explicit. An important example of this family of methods is the convex-concaveadditive decomposition of the free energy introduced by Eyre [50] for general gradient flows, particularly for the CHequation. In this case, since the gradient term is quadratic, it contributes to the convex part of the decomposition, thechallenge being the decomposition of Ψ , the bulk free energy, generally a non-convex function like a double-well.The additive decomposition of Ψ is done by splitting concave and convex terms of the functional such that the convexterms are treated implicitly, and the concave terms can remain explicit. However, the additive decomposition of Ψ isnot unique for all the possible functions of Ψ [49]. The splitting method proposed in [50] yields an unconditionallyenergy-stable method and is unconditionally uniquely solvable, although the proved local truncation error is justsecond-order accurate in the timestep size, rendering a first-order accurate in time method. Some studies developedsecond-order convex-splitting time integration schemes for the CH equation [51, 32]. However, unconditionally-energystability and unconditionally uniquely solvability properties are not yet achieved for a general form of Ψ , within the freeenergy decay original approach or without the use of numerical stabilization [49]. The development of second-order,unconditionally-energy stable, and unconditionally uniquely solvable integration methods is still an active researchtopic in the present context.Vignal et al. [33] introduced another approach to derive second-order unconditionally energy-stable time integrators forphase-field models. The main idea is to use Taylor’s series expansions of the bulk free energy function Ψ to derivean approximation for Ψ (cid:48) = ∂ Ψ ∂φ . This method is mathematically proven unconditionally energy-stable regardless ofmesh and time step size and second-order accurate in time for quartic potentials. Although there is no proof that thismethod is unconditionally uniquely solvable, it works well with adaptive time-stepping [33]. Due to its simplicity and4 PREPRINT - O
CTOBER
1, 2020desirable properties, this is the method chosen in this study. Applying Vignal et al. [33] time integration method tothe semi-discrete variational formulation of the problem given in equation (9), where the subindex n is the time stepnumber and considering the initial conditions φ ( x ,
0) = φ , the fully discrete system is as follows: (cid:18) w h , [[ φ ]]∆ t n +1 (cid:19) Ω + ( ∇ w h , M ( φ h ) ∇ µ hn +1 ) Ω + ( w h , σ { φ } ) Ω − ( w h , σ ¯ φ h ) Ω = 0 , (10) ( q h , µ hn +1 ) Ω − ( q h , ˜Ψ (cid:48) ) Ω − ( ∇ q h , (cid:15) ∇{ φ } ) Ω = 0 , (11)where [[ φ ]] = φ hn +1 − φ hn , ∆ t n +1 = t n +1 − t n , { φ } = φ hn +1 + φ hn , and the approximation ˜Ψ (cid:48) is defined as ˜Ψ (cid:48) = ∂ Ψ n +1 ∂φ − ∂ Ψ n +1 ∂φ [[ φ ]]2 + ∂ Ψ n +1 ∂φ [[ φ ]] . (12)By being fully implicit, this scheme suits the proposed adaptive time-stepping strategy, allowing larger time stepswithout compromising the stability of the method. It is essential to mention that the Taylor expansion of Ψ presented inEq. (12) has guaranteed convergence for quartic potentials representing the homogeneous free energy function. In thisstudy, we apply the time integration scheme for the NCH equation derived from the O–K functional. Despite beingmathematically designed for the CH equation, the employment of the described temporal integration method on theNCH equation holds the properties of free energy decay and mass conservation, as seen in Section 5. The choice of a proper time integration method for the NCH equations is a difficult task, since, in several physicalsituations, the equation has different time scales, creating a tradeoff between accuracy and performance. For instance,the initial stage of the phase segregation of a mixture is dictated by fast dynamics, requiring small step sizes while thelatter stages reveal slow dynamics, allowing large time steps. Thus, to improve the efficiency of the computations, atime adaptivity scheme is often used to automatically change the time step size to capture both fast and slow dynamicsof the equations, as well as the nonlocal dynamics inherent to the NCH equation, improving the performance of thesimulations without any accuracy loss.Studies in the literature discuss time adaptivity schemes for the CH equation. Some schemes rely on the evaluation ofthe Ginzburg-Landau free energy [52, 53], requiring, however, the tuning of very sensitive empirical parameters [53].Different approaches are seen in [32, 34, 31, 33, 51, 35] are based on simple time-step controllers and present morerobustness and better results. In this study, we assess different controllers applied to the NCH equation.
The NCH equation can be expressed as a dynamical system of the form: ˙ φ = F ( φ ) ,φ (0) = φ , (13)where φ ∈ R n sd and F : R n sd → R n sd is a Lipschitz map. Since the time integration method used in this study is aone-step method, considering a step size ∆ t , there is a map Φ : R n sd → R n sd such that: φ n +1 = Φ( φ n ) ,φ (0) = φ . (14)Equation (14) is a discrete-time dynamical system that approximates Equation (13). It is possible to use the sameapproach for an additional map Ξ : R → R to vary the step size: ∆ t n +1 = Ξ(∆ t n ) . (15)where ∆ t n +1 = t n +1 − t n and ∆ t n = t n − t n − .The map Ξ uses information about the numerical solution φ n when defining the new step size ( ∆ t n +1 ) while the map Φ is based on the time step ∆ t n +1 . An adaptive time–stepping method can be expressed as the following recursions: φ n +1 = Φ( φ n ) , ∆ t n +1 = Ξ(∆ t n ) . (16)5 PREPRINT - O
CTOBER
1, 2020Figure 1: Adaptive time-stepping viewed as a feedback control system. Adapted from [55].We assume that the relation between the error and the step size is asymptotic, that is: r n = | ζ n | ∆ t κn (17)where r n is the norm of the local error estimate, | ζ n | is the norm of the principal error function, and κ is related to theorder of the method. In our case, the principal error function can be viewed as a disturbance in the system, such as aNewton solver residual, and the integration method is second-order accurate, so κ = 2 . It can be seen that r n → if ∆ t → .The idea behind the use of control theory on adaptive time-stepping is that the map Ξ controls an estimated numericalerror within a prescribed tolerance, T OL . The mathematical background containing the detailed description of theuse of control theory in temporal adaptivity in ODEs is given in [54, 55, 56]. The recursion can be translated into aclosed-loop, a common dynamic structure in Control Theory, as seen in Figure 1. In this sense, the use of a linearfeedback controller on step size adaptivity is translated in: given a present time step t n , the controller defines a futuretime step ∆ t n +1 such that the local error of the present time step r n is controlled within a given tolerance T OL by acontroller whose properties and tuning are defined in the mapping Ξ . The time step is evaluated within the process(which is solving the NCH equation in that given time) and feeding the error r n +1 to the controller, restarting theloop. If the estimated error is not within the prescribed tolerance, the controller reevaluates a new, smaller time stepsize until this condition is satisfied. Besides the expectancy of reducing the number of linear/nonlinear systems to besolved, Control Theory provides smoother step size sequences (which improves the solution regarding smoothness[54]), improved computational stability, and a regular, tight tolerance proportionality. The literature presents various methods to estimate the temporal error on ordinary differential equations (ODEs). Theinterested reader can see in [57, 58]. The usual strategy to estimate the temporal error is to evaluate an ODE solution ata given time step with integration methods of a different order of accuracy and compute the norm of the difference ofthe solutions relative to the norm of the solution obtained by the higher-order method. This strategy is seen in the CHcontext in [32, 34, 31] and it can be mathematically written as, r = || φ n +1 − ˆ φ n +1 |||| φ n +1 || , (18)where φ n +1 is obtained through an integration method of a higher order than the solution obtained in ˆ φ n +1 . Althoughthis is considered a common error estimator, the calculation of the time step n + 1 twice is required to obtain φ n +1 and ˆ φ n +1 . In [35], φ n +1 is obtained by taking two time steps of size ∆ t n +1 / while ˆ φ n +1 is obtained by considering thetime step ∆ t . In this case, the authors use the same integration method for both solutions, but each time step needs tobe computed three times.A different approach is seen in [33]. Consider an error estimation where the solutions at t n and t n − are stored andthe error is estimated a posteriori by extrapolation by a lower-order time integration method, since the solution t n +1 is obtained with a second-order scheme. This estimation is done through variable step-size backward differentiation,where the error obtained in the lower-order method is controlled. In the present work, we consider the lower-ordermethod to be the backward-Euler method. Therefore, the local truncation error of the backward-Euler method is: τ BE ( t n +1 ) = − ∆ t φ (cid:48)(cid:48) ( t n +1 ) + O (∆ t ) . (19)Given the stored solutions φ n +1 , φ n and φ n − at times t n +1 , t n and t n − respectively and neglecting the effects of the O (∆ t ) terms, equation (19) can be approximated by the variable step-size backward difference formula. So the error6 PREPRINT - O
CTOBER
1, 2020Table 1: Controllers parameters.Controller κ P κ I κ D κ T Integral (I) . . . . PID .
075 0 .
175 0 .
01 0 . PC11 .
333 0 .
333 0 . . estimation is now: E n +1 = − η φ n +1 + 1 η − φ n − η ( η − φ n − (20)where η = (∆ t n +1 + ∆ t n ) / ∆ t n +1 .With the error function, the weighted local truncation error (WLTE) can be written as: r = (cid:118)(cid:117)(cid:117)(cid:116) n nodes n nodes (cid:88) i =1 (cid:18) E in +1 τ abs + τ rel max ( | φ in +1 | , | φ in +1 + E i ( n +1) | ) (cid:19) (21)where τ abs and τ rel define tunable absolute and relative tolerances, respectively, and the index i = 1 , , ...n nodes refersto the nodal index.The weighted local truncation error r is used to control the error at each time step. By definition, values of r ≤ mean that the local truncation error is within the user-prescribed tolerances. In this case, the step just taken canbe accepted, and the time integration can move forward with either the same or a larger time step size. On thecontrary, values of WLTE larger than one imply unacceptable errors. That said, the step taken is then rejected andretaken with a smaller time step size. We also constrain the time step such that ∆ t n +1 ∈ [∆ t min , ∆ t max ] , thatis, ∆ t n +1 = max (∆ t n +1 , ∆ t min ) and ∆ t n +1 = min (∆ t n +1 , ∆ t max ) , where ∆ t min and ∆ t max are user suppliedparameters. There are several timestep controllers in the literature and many ways to classify them [55]. In this study, we considerthree controllers: an integral controller, a PID controller, and the PC11 predictive controller. The three controllers haveunique properties and have been used in the context of time adaptivity of PDEs such as the Navier-Stokes [59] and theconvection-diffusion equation [60]. The integral controller is the simplest and controls the relationship between theerror in the present and past time. This simplicity is known to grant the integral controller a large number of rejectedsteps [31]. The PID controller has three controlling terms - proportional, integral, and derivative - that adjust the timestep to changes in the estimated errors of the last three time steps. The predictive controller PC11 is suggested for timeadaptivity in stiff equations [54], which is the case for the CH/NCH equations and has a different structure compared tothe other two controllers. The three controllers can be written as: ∆ t n +1 = ρ (cid:18) r n r n +1 (cid:19) κ P (cid:18) r n +1 (cid:19) κ I (cid:18) r n r n +1 r n − (cid:19) κ D (cid:18) ∆ t n ∆ t n − (cid:19) κ T ∆ t n (22)where the parameter ρ is a safety factor used to smooth the time step growth. In the literature, it is common to see ρ = 0 . , although in [51] other values for ρ and tolerances for other phase-field computations were proven better [51].Here, we adopt ρ = 0 . unless stated otherwise. We evaluate these parameters for a nonlocal case in the numericalvalidation section and investigate the accuracy and performance results. The parameters κ P , κ I , κ D and κ T for eachcontroller are shown in Table 1. To avoid tuning the controller parameters, which can be very time consuming, theparameters for the I controller [54], the PID controller [59] and the PC11 controller [60] are taken from the literature. Remark:
Although the use of time adaptivity schemes based on the linear feedback control theory has not yet beenexplicitly mentioned in the CH equation literature, the integral controller has been used by others [31, 32, 33, 34, 51].Even when the PID controller is used for CH equation [31], the integration method for the PID error estimation is notguaranteed to be energy stable. The use of a predictive controller for the NCH equation is unprecedented. Moreover,the error estimation is based on solution norms of time integration methods with different accuracy, requiring thecalculation of the same step twice and, therefore, can be time-consuming. In the present work, we employ an errorestimation method based on extrapolation that avoids computing the same time step twice.7
PREPRINT - O
CTOBER
1, 2020Figure 2: Time history of the time step ∆ t during the 2D spinodal decomposition simulation. Numerical simulations are made to validate the time adaptivity strategy. Initially, we consider a case with no-fluxboundary conditions, ¯ φ = 0 . and σ = 0 , that is, a standard CH simulation of spinodal decomposition. We considerfour simulations: one for each controller, and a fixed time step simulation. We compare the frames in all simulationsto check if they all represent the same physical stages. We consider a square domain with a nodes. The squaredomain is divided into cells, each cell discretized by two linear triangles. For our simulations, we consider τ abs = τ rel = 10 − , Ψ( φ ) = ( φ − , M ( φ ) = 1 and the initial time step ∆ t = 1 × − . We also constrain thetime step size to the limits ∆ t min = 1 × − and ∆ t max = 5 × − . Preliminary simulations showed that whenthe simulation approaches the steady-state, the controller allows time step sizes of O (10 − ) . Although accepted by thecontrollers, the use of time step sizes of this magnitude increases the number of linear iterations significantly in the laterstages, increasing the computational cost. We also note that the number of rejected steps for the three controllers aresignificantly reduced, suggesting that limiting the time step improves the controllers’ behavior. The FEniCS frameworkv2019.1.0 invokes several linear algebra backend packages to solve the linear and nonlinear systems arising from thefinite element method. In our case, we choose the PETSc package, inheriting Newton’s method to solve the nonlinearsystems with a relative tolerance η NL = 10 − and the GMRES solver with Block-Jacobi ILU(0) preconditioner withrelative tolerance η r = 10 − and absolute tolerance η a = 10 − for the linear systems.Figure 2 shows the time step history for the four simulations. We can see in this figure the multiscale nature of the CHequation in the spinodal decomposition. In the initial stage where the phases are being defined, The time step increasesby orders of magnitude and then decays (from t = 10 − to t = 10 − , approximately). This strong variation occurs dueto the rapid dynamics in the early stages of the spinodal decomposition. Consequently, the controller produces smallertime steps to keep the estimated error within the prescribed tolerance. The intermediate stage begins when all bubbles8 PREPRINT - O
CTOBER
1, 2020Figure 3: Free energy for the three simulations and some snapshots describing phase separation.are approximate of the same size, leading to simultaneous Ostwald’s ripening events in the domain, preventing thetime step from growing. The controller keeps the oscillations in this interval (from t = 10 − to t = 10 − ) to capturethe bubble shrinkage. The final stages are where the time step has larger values, indicating slow dynamics involvingsurface motion. The final stage, however, still requires small time steps when an Ostwald ripening is occurring, but thecontroller allows the time step size to grow by several orders of magnitude since no rapid dynamics are seen in thisstage when there is no shrinkage. A few simulation snapshots at the different stages and the free energy decay for thethree controllers are seen in Fig. 3.We now evaluate the use of time step controllers in terms of physical accuracy and performance. In terms of physicalaccuracy, we observe from both Figs 2 and 3 that the three curves regarding the adaptive time stepping simulationsare practically overlapping in terms of free energy and time step size. The overlapping is a strong indicator that thesimulations are practically identical, meaning that the results are physically consistent for all three controllers. Figure 4presents the results for all the simulations at t = 0 . ± . , showing that all the controllers yield the samesolution.We note that the early stages of the spinodal decomposition demand a smaller time step size than the maximum allowedfor the fixed time step simulation. This fact means that the early stages for the fixed time step simulation are obtainedwith a smaller amount of time steps than the adaptive simulations. However, this is valid only to the initial simulationstage. Afterward, adaptive simulations become more efficient. Figure 5 shows how many time steps are needed foreach method to reach the point where the fixed time step is no longer more efficient. We observe that the PID demandsmore time steps to reach the same simulation stage as the fixed time step simulation compared to the integral and PC11controllers. It is expected that the PID controller behaves more conservatively, since its formulation carries the estimatederror in three different time steps, making it a more rigid controller than the others. Table 2 show a comparison of9 PREPRINT - O
CTOBER
1, 2020 (a) Fixed time step size. (b) Integral controller.(c) PID controller. (d) PC11 controller.
Figure 4: Comparison of the four simulations at t = 0 . ± . .10 PREPRINT - O
CTOBER
1, 2020Figure 5: Time steps required for each controller to reach the same physical stage as the fixed time step simulation.the performance results for the three controllers, and Figure 6 presents the evolution of linear and nonlinear iterationsduring the simulation.All methods significantly improved the spinodal decomposition simulation since it is possible to reach larger simulationtimes with a smaller number of steps. We observe from Fig. 3 that the steady-state is reached at around t = 0 . in oursimulations. To reach the steady-state using the fixed time step scheme, considering that ∆ t = 2 × − is the largestpossible fixed time step that would not introduce unacceptable errors in the simulation, it would be necessary . × time steps, while using the time step adaptivity, it is reached with circa , time steps, as seen in Table 2. In relativeterms, adaptive simulations reach the steady-state in approximately . of the simulation time needed for a fixed timestep simulation, reinforcing the importance of temporal adaptivity. Observe that, in terms of required time steps, thePC11 required fewer time steps than the other two controllers while the PID controller solved the larger amount of timesteps. The Integral controller, the simplest controller, has less control over the growth of the time steps, presentingmore rejected steps. We also evaluate the performance in terms of linear and nonlinear iterations. We consider theabsolute CPU effort calculated as the total number of linear iterations during the simulation, considering accepted andrejected steps. The controller with a larger absolute CPU effort becomes the reference for calculating the relative CPUeffort. Comparing in Table 2 and Figure 6 the three adaptive simulations, the PC11 controller has the smaller number oftotal linear iterations with an improvement of over the amount of the same quantity for the PID controller and in comparison with the I controller. Even though the PID solution presents the lower average of nonlinear andlinear iterations, it requires more time steps. By computing the total CPU effort, that is, the total number of iterationsevaluated, we see that the other two controllers have a better performance.11 PREPRINT - O
CTOBER
1, 2020Table 2: Performance results for the time adaptivity schemes for each time step controller in the 2D spinodal decompo-sition. Step size Accepted Rejected Avg. Nonlinear Avg. linear Relative CPUController Steps Steps Iterations Iterations EffortI . . . PID . . . PC11 . . . Integral PID PC11 Integral PID PC11
Figure 6: Number of nonlinear (left) and linear iterations (right) for each time step during the adaptive simulations.Rejected steps included. Solver tolerances: η NL = 10 − , η r = 10 − , η a = 10 − . In this section, we solve the NCH equation in two and three dimensions to evaluate the performance and accuracyof the time step controllers. The same parameters regarding the domain, mesh size, interface thickness, free energyhomogeneous function, and other numerical parameters are extended from the previous examples to the nonlocal cases.The boundary conditions, however, are considered periodical to preserve the NCH equation pattern formation. It isknown that the variation of the parameters (cid:15) , σ , ¯ φ , and the domain size interfere directly with the steady-state structureof the NCH equation. We define these parameters such that the minimizers are situated on a locally stable region of thephase diagram and a domain size large enough compared to the intrinsic length scale of the minimizers of the O-Kfunctional [5]. Initially, we consider a 2D case, where the copolymers steady-state has a hexagonally packed spotsstructure, as seen in the phase diagrams in [27, 28]. We consider ¯ φ = 0 . , (cid:15) = 0 . and σ = 500 . For this first nonlocalexample, we consider an assessment of the controllers’ parameters. In [51], there is a remark that the use of controllersfor time step size adaptivity for the Swift-Hohenberg equation [61, 62] with a smaller safety coefficient and tightertolerances yields a smaller percentage of rejected time steps and, according to [28], the NCH can be viewed as a hybridof the Swift–Hohenberg equation and the CH equation. Therefore, we compare the results of the standard controllerparameters with the results obtained by considering ρ = 0 . and τ abs = τ rel = 1 × − . We label the simulationsfor the controller parameters used in the CH example as Case 1 and Case 2 for the new proposed values. Figure 7 showsthe time step size history for the six simulations, while the corresponding steady-state configurations are seen in Figure8. Table 3 shows the results for the controllers’ performance in all cases.We note in both cases in Fig. 7 that the time step size curves are not overlapping, as in the previous examples. Therefore,for a given instant where the curves do not match, the phenomenon in Fig. 4 is not observed in the diblock copolymercontext, meaning that different controllers lead to different observations in the time evolution of the diblock copolymermelt. However, this does not affect the formation of the steady-state structures related to the selected set of parameters,which is the information of interest in most cases. We observe the steady-state for all six simulations of Fig. 8 and notethat the six simulations converged to hexagonally packed spots, as initially predicted by the phase diagrams [27, 28],despite minor differences that arise due to the periodic boundary conditions. We also observe in Figure 9 that oursimulations do not present any unphysical properties in the free energy decay and the mass conservation for the bestperforming controllers. 12 PREPRINT - O
CTOBER
1, 2020
Integral PID PC11 FixedIntegral PID PC11 Fixed
Figure 7: Time step size history for Case 1 (top) and Case 2 (bottom).13
PREPRINT - O
CTOBER
1, 2020 (a) I controller. (b) I controller.(c) PID controller. (d) PID controller.(e) PC11 controller. (f) PC11 controller.
Figure 8: Steady-state for the nonlocal Cahn-Hilliard equation in two dimensions for the three controllers for Case 1(left) and Case 2 (right). 14
PREPRINT - O
CTOBER
1, 2020Table 3: Results for the time adaptivity schemes for each time step controller in the NCH simulations.Time Step Accepted Rejected Avg. Nonlinear Avg. Linear Relative CPUController Steps Steps Iterations Iterations EffortI . . . Case 1 PID . . . PC11 . . . I . . . Case 2 PID . . . PC11 . . . PC11 PC11
Integral Integral
Figure 9: Free energy decay and mass conservation for the best performing controllers for Case 1 (top) and Case 2(bottom) simulations.Regarding the performance, we can see in Table 3 that reducing the prescribed tolerances leads to a significant increasein the required number of time steps to reach the steady-state while the rejected steps decreased. We observe thatthe best performing controller in our simulations is PC11 for Case 1, while in Case 2, the I controller shows betterperformance. One possible explanation is that the I controller’s aggressive behavior combined with the tighter tolerancesin Case 2 leads to a more controlled environment where the number of rejected steps is not significant compared to theincreased number of time steps required for the simulations to reach the steady-state. Nevertheless, for both cases, the Icontroller presents the largest average number of linear iterations while the PID controller has the smaller.After observing the effects of the controllers’ parameters on the NCH equation, we note that the use of the standardvalues used in Case 1 requires less computational effort and does not influence the steady-state evaluation. Therefore,we extend our analysis using Case 1 for different parameter sets. We consider three examples: test case A, where ¯ φ = 0 . and σ = 1000 , test case B where ¯ φ = 0 . and σ = 500 and test case C where ¯ φ = 0 . and σ = 1000 . Figure15 PREPRINT - O
CTOBER
1, 2020 (a) Test case A (b) Test case B(c) Test case C
Figure 10: Steady-state for the NCH equation in two dimensions.10 shows the steady-state these test cases. The minimizing structure of the melt in two dimensions are hexagonallypacked spots, stripes and mixed states [28, 27]. In the figure, we can see spots (Fig. 10a), stripes (Figs. 10b) and mixedstructure (Figs. 10c) melts. All the simulations reached the steady-state at around t = 0 . .Figure 11 shows the time step size history for these simulations. We see that in test case A the PC11 controller presentsa sharp decrease in the time step size in the last simulation stages, but as the simulation approaches the steady-state,the time step size increases, returning to values of the same order of the other two controllers. Table 4 and Figure 12show for all test cases the performance results for the three controllers. The I controller is the best in test cases A and B,while PC11 is the best performing controller in test case C. We also notice a smaller number of accepted and rejectedtime steps for simulations with ¯ φ = 0 . , that is, test cases B and C, in comparison with simulations where ¯ φ = 0 . (test case A). The number of linear iterations is smaller for the three controllers in case A, but the number of nonlineariterations is of the same order in all test cases for the three controllers. The PC11 controller has the larger averagenumber of linear iterations on Test Case A, but the same behavior from the previous simulations is seen on Test Cases Band C, which is, the I controller has the largest number of average linear iterations while the PID presents the least.We note in Figure 12 that the solutions with PC11 present a sharp increase in the number of linear iterations around t = 1 × − , particularly for test case A. The PID solutions at the same time interval exhibit the lower number ofiterations. In all three test cases, the number of linear iterations increases when we approach the steady-state. We then16 PREPRINT - O
CTOBER
1, 2020
Integral PID PC11 Fixed (a) Test case A
Integral PID PC11 Fixed (b) Test case B
Integral PID PC11 Fixed (c) Test case C
Figure 11: Time step size history for the cases described in Fig. 10.17
PREPRINT - O
CTOBER
1, 2020
Integral PID PC11 Integral PID PC11Integral PID PC11 Integral PID PC11Integral PID PC11 Integral PID PC11
Figure 12: Number of nonlinear (left) and linear iterations (right) for test cases A, B and C, respectively. Rejected stepsincluded. Solver tolerances: η NL = 10 − , η r = 10 − , η a = 10 − .18 PREPRINT - O
CTOBER
1, 2020Table 4: Results for the time adaptivity schemes for each time step controller in the 2D diblock copolymer simulations.Test Time Step Accepted Rejected Avg. Nonlinear Avg. Linear Relative CPUCase Controller Steps Steps Iterations Iterations EffortI . . . A PID . . . PC11 . . . I . . . B PID . . . PC11 . . . I . . . C PID . . . PC11 . . . Table 5: Results for the time adaptivity schemes for each time step controller in the 3D diblock copolymer simulations. ¯ φ Time Step Accepted Rejected Avg. Nonlinear Avg. Linear Relative CPUController Steps Steps Iterations Iterations EffortI . . . . PID . . . PC11 . . . I . . . . PID . . . PC11 . . . evaluate the free energy decay and mass conservation, as shown in Fig. 13 for the best performing controllers. We seethat the free energy decays for all test cases. Mass is for practical purposes conserved, with losses of order − . Recallthat the tolerance for the nonlinear solver is η NL = 10 − , and for the linear solver η r = 10 − and η a = 10 − . Thusthe values for mass conservation are compatible with the accuracy obtained in each time step solution.We extend our analysis to three dimensions. We consider a cubic domain with nodes, trilinear hexahedral elements, σ = 500 and evaluate the cases where ¯ φ = 0 . and . . Figure 14 shows the melt structure for the two 3D simulationswith the consistent patterns observed in this phenomenon. The chosen parameters lead to stable melts, where the casewhere ¯ φ = 0 . is a bicontinuous melt and ¯ φ = 0 . leads to a perforated layer melt [63]. Steady-state is reached onapproximately t = 1 . for the case where ¯ φ = 0 . and t = 0 . for ¯ φ = 0 . . Initially, we discuss the time step sizehistories shown in Figure 15. We can see that for ¯ φ = 0 . and ¯ φ = 0 . , the time step histories exhibit an initial stagewith a fast time step-growth, an intermediate stage, where the time step increases but oscillates, and the final stage wherethe time step recover a fast growth. For ¯ φ = 0 . we also see smaller oscillations in the final stage. We note that thecase where ¯ φ = 0 . is different from the 2D test cases B and C. In the 3D case, the time step size reveals a much morecomplex behavior than the 2D case. This behavior is related to the generation of the complex structure melts exhibitedin Figure 14. In terms of efficiency, Fig. 16 and Table 5 show the performance data and the time history of the number ofnonlinear and linear iterations for the 3D cases. We note that the number of accepted and rejected steps and the averagenumber of nonlinear and linear iterations increased compared to the 2D test cases, reflecting the higher complexityexistent in 3D copolymer simulations. We observe that the I controller has the largest number of linear iterations againwhile the PID controller the smaller. Also, we can see in Fig. 16 that for both cases ( ¯ φ = 0 . and ¯ φ = 0 . ) the numberof linear iterations exhibits a remarkable growth as the solution approaches the steady-state. This growth is exceedinglylarge for the I and PID controllers. We observe that the PC11 controller exhibits the best performance, saving around of the computational effort required by the I and PID controllers for the case where ¯ φ = 0 . and almost a gain for ¯ φ = 0 . . In terms of accuracy, we observe the free energy decay and mass conservation properties in Fig. 17for the best performing controllers for each case. Again, mass conservation is within the accuracy obtained at each timestep solve. This paper presents time adaptivity schemes for the nonlocal Cahn-Hilliard equation derived from the Ohta-Kawasakifree energy functional. We have unified the time adaptivity schemes under the linear feedback control theory. Theerror estimate for the time adaptivity schemes results from an extrapolation method proposed in [33], avoiding solvingthe problem twice at the same time step. We test our scheme on simple examples such as a 2D phase separation with19
PREPRINT - O
CTOBER
1, 2020
Integral IntegralIntegral Integral
PC11 PC11
Figure 13: Free energy and mass conservation for the test cases A, B and C, from top to bottom, for the best performingparameters. 20
PREPRINT - O
CTOBER
1, 2020Figure 14: Structure of a monomer on the generated diblock copolymer melts from the 3D simulations for ¯ φ = 0 . (top) and . (bottom). 21 PREPRINT - O
CTOBER
1, 2020
Integral PID PC11 FixedIntegral PID PC11 Fixed
Figure 15: Time step size history for the three controllers. Top, ¯ φ = 0 . and bottom ¯ φ = 0 . PREPRINT - O
CTOBER
1, 2020
Integral PID PC11 Integral PID PC11Integral PID PC11 Integral PID PC11
Figure 16: Number of nonlinear and linear iterations during simulation for the three controllers. Rejected steps included.Top ¯ φ = 0 . ; bottom, ¯ φ = 0 . . Solver tolerances: η NL = 10 − , η r = 10 − , and η a = 10 − .constant mobility and more challenging ones, the diblock copolymer self-assembly in two and three dimensions. Weevaluate the accuracy and performance of three time step size controllers. The simulations reveal the PID controller’sconservative behavior, which required the computation of more times steps than the other controllers to reach thesteady-state and the I controller’s aggressive behavior, with many rejected time steps. Our simulations for the CHequation suggest that the PC11 controller is the best performing controller. However, for the 2D NCH simulations, theresults vary. Initially, we assessed the values for the tolerance and safety coefficient and observed that the I controllerpresents the best performance compared to the other controllers when the tolerance is stricter and the safety coefficientsmaller. For the standard values used in the literature, the results obtained are accurate, and the PC11 yields the bestperformance. Besides, we consider different copolymer parameters and evaluate their influence on the controllers’performance. We note that in some cases, the PC11 and the I controllers have a better performance. Subsequently, weconsider two 3D NCH simulations. In both, PC11 is the most efficient controller. Furthermore, in terms of accuracy,our numerical results with the second-order energy stable time integration method introduced in [33] coupled withtemporal adaptivity show numerical of evidence mass conservation and free energy decay for the nonlocal case, asanticipated theoretically in [37].We note that the controllers reveal a subtle interplay between the time step size and the nonlinear and linear systemsolves. For most of the simulations presented in this study, the I controller has the largest average number of lineariterations while the PID controller the least. This effect is an important metric in the sense that the average numberof linear and nonlinear iterations must be taken into account with the total number of times steps to evaluate theperformance of a time step size controller. For instance, for the 3D case where ¯ φ = 0 . , the I controller presents asmaller number of total time steps (accepted and rejected) than the PC11, but the latter has a smaller overall numberof linear iterations. This observation is important because we observe a performance gain by using the PC11controller in this time and resource-demanding simulation. In terms of accuracy, we also notice that the time adaptivityscheme reproduces all different physics, such as those found in the local and the nonlocal Cahn-Hilliard equation in twoand three dimensions. An essential statement obtained from this study is that one should always consider temporaladaptivity for the nonlocal Cahn-Hilliard equations since the gains in the computational effort are substantial.23 PREPRINT - O
CTOBER
1, 2020
PC11 PC11PC11 PC11
Figure 17: Free energy and mass conservation for the 3D simulations for the best performing controllers. Top ¯ φ = 0 . ;bottom, ¯ φ = 0 . . Acknowledgements
This research was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil(CAPES) - Finance Code 001. This research has also received funding from CNPq and FAPERJ. Computer time inLobo Carneiro supercomputer was provided by the High Performance Computer Center at COPPE/Federal Universityof Rio de Janeiro, Brazil.
References [1] Cahn JW, Hilliard JE. Free energy of a nonuniform system. I. Interfacial Free Energy.
The Journal of ChemicalPhysics
The Journal of ChemicalPhysics
Mathematical Problems in Engineering
Macromolecules
SIAM Journal on Applied Mathematics
IEEETransactions on Image Processing
PREPRINT - O
CTOBER
1, 2020[7] Hohenberg P, Halperin B. Theory of dynamic critical phenomena.
Reviews of Modern Physics
Mathematical Models and Methods in Applied Sciences
Proceedingsof the Royal Society A: Mathematical, Physical and Engineering Sciences
Mathematical Models and Methods in Applied Sciences
Journal of Fluid Mechanics
Computer Methods in Applied Mechanics and Engineering
Journal of theMechanics and Physics of Solids
Mathematical and Computer Modelling
International Journal for Numerical Methodsin Biomedical Engineering
Structural and Multidisciplinary Optimization
Journal of Computational Physics
Journal of Statistical Physics
SIAM Journal on Applied Mathematics
Journal of Mathematical Analysis and Applica-tions
The Physics of Block Copolymers . 19. Oxford University Press: New York, NY, USA . 1998.[22] Kim H, Park S, Hinsberg W. Block Copolymer Based Nanostructures: Materials, Processes, and Applications toElectronics.
Chem. Rev.
Journal of the American Mathematical Society
Interfaces and Free Boundaries
Chaos
Nonlinear Analysis
Nonlin-earity
SIAM Journal on Applied Dynamical Systems
Current Applied Physics
SIAM Journal on Matrix Analysis andApplications
Journal of Computational Physics
PREPRINT - O
CTOBER
1, 2020[32] Gómez H, Calo VM, Bazilevs Y, Hughes TJR. Isogeometric analysis of the Cahn-Hilliard phase-field model.
Computer Methods in Applied Mechanics and Engineering
Computer Methods in Applied Mechanics and Engineering
Journal of ComputationalPhysics
International Journal for Numerical Methods inEngineering
Numerical Approximation of the Ohta-Kawasaki Functional . M.Sc. thesis. Kellogg College, Universityof Oxford, UK; 2012.[37] Gal CG. Doubly nonlocal Cahn–Hilliard equations.
Annales de l’Institut Henri Poincare (C) Analyse Non Lineaire
Applied Mathematics Letters
CRM Proceedings and Lecture Notes: Singularities in PDE and the Calculus ofVariations, American Math. Society
Archive for Rational Mechanics and Analysis
MRS Proceedings
SIAM Journal onNumerical Analysis
NumerischeMathematik
Applied NumericalMathematics
Journal ofComputational Physics
Numerical Algorithms
Automated Solution of Differential Equations by the Finite Element Method .Springer . 2012.[48] Alnæs MS, Blechta J, Hake J, et al. The FEniCS Project Version 1.5.
Archive of Numerical Software
Computational Phase-Field Modeling
Unpublished article
Applied Numerical Mathematics
Computers and Mathematics with Applications
Communications in Computa-tional Physics
Numerical Algorithms
ACM Transactions on Mathematical Software
Applied NumericalMathematics
PREPRINT - O
CTOBER
1, 2020[57] Hairer E, Nørsett S, Wanner G.
Solving Ordinary Differential Equations I Nonstiff Problems . Berlin: Springer-Verlag Berlin Heidelberg . 1993.[58] Hairer E, Wanner G.
Solving Ordinary Differential Equations II Stiff and Differential-Algebraic Problems . Berlin:Springer-Verlag Berlin Heidelberg . 1996.[59] Valli AMP, Carey GF, Coutinho ALGA. Control strategies for timestep selection in finite element simulation ofincompressible flows and coupled reaction-convection-diffusion processes.
International Journal for NumericalMethods in Fluids
Computer Methods in Applied Mechanics and Engineering
Reviews of Modern Physics
Physical Review A