{'fraction_non_alphanumeric': 0.04433630711002974, 'fraction_numerical': 0.016761286834279535, 'mean_word_length': 4.682795698924731, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Improved time integration for phase-field crystal models of solidification

We optimize a numerical time-stabilization routine for the phase-field crystal (PFC) models of solidification. By numerical experiments, we showcase that our approach can improve the accuracy of underlying time integration schemes by a few orders of magnitude. We investigate different time integration schemes. Moreover, as a prototypical example for applications, we extend our numerical approach to a PFC model of solidification with an explicit temperature coupling.


Introduction
The so-called phase-field crystal (PFC) model [1][2][3] is a prominent approach to describe crystal structures at large (diffusive) timescales through a continuous, periodic order parameter representing the atomic density. The PFC model is based on a Swift-Hohenberg free energy functional [1][2][3], which can be written as The scalar order parameter ψ(r, t) 1 is related to the atomic number density, and Ω ∈ R n (n = 2 in this work). A set of parameters λ, κ > 0 characterize the phase space and material properties. Together with appropriate boundary-and initial conditions, the dynamics of ψ is described via a conservative (H −1 ) gradient flow of F . The time evolution reads with L[ψ] and N [ψ] a linear and nonlinear differential operator, respectively: Numerical solutions of the partial differential equation (2) can be computed using a Fourier pseudo-spectral method for spatial discretization, enforcing periodic boundary conditions in combination with a suitable time discretization [8][9][10][11]. As time stepping schemes, we consider a linear first-order semi-implicit (IMEX) scheme, which treats L[ψ] fully implicitly and N [ψ] fully explicitly, a first-order exponential time integration method (ETD1), and a second-order exponential time integration method of Runge Kutta type (ETD2RK) [4].
Additionally, we exploit the numerical time-stabilization routine presented in [5]. In combination with an underlying time stepping scheme, this approach features a convex-concave splitting of the free energy controlled by parameters C 1 and C 2 : We show in the following that they can be tuned to achieve an improved numerical accuracy, and that the best performances are here obtained for negative values. For every time integration scheme, we determine optimal C 1 , C 2 by minimizing the difference of the error (least square of the difference) in the free energy decay with respect to a numerical reference solution and denote the resulting time stepping schemes by C-IMEX, C-ETD1, and C-ETD2RK.

Numerical parameter study
A convergence study of the considered numerical methods is reported in Figure 1, for a specific set of model parameters leading to the growth of a crystal in a domain Ω = [−200, 200] 2 (initial condition as in Figure 2 (a)). The numerical simulations are performed on an uniform grid with an element size of dx = dy = 0.78125 and for time steps ∆t ∈ [10 −2 , 5]. Figure 1 shows the residual of ψ, R ψ , evaluated as the discrete L 2 distance from a numerical reference solution for different ∆t.
The schemes IMEX, ETD1 and ETD2RK (C 1 = C 2 = 0) and the corresponding ones including the splitting (4), C-IMEX, C-ETD1, and C-ETD2RK are compared, with a reference solution corresponding to the C-ETD2RK scheme with ∆t = 0.01. As expected, the (C-)IMEX and (C-)ETD1 schemes converge linearly, and the (C-)ETD2RK scheme converges quadratically for a decreasing time step size ∆t. For a fixed ∆t, the C-IMEX and C-ETD1 schemes give two orders of magnitude smaller errors than the IMEX and ETD1 schemes, respectively, whereas the C-ETD2RK scheme gives a half order of magnitude smaller error than the ETD2RK scheme. By fixing an accuracy instead, the C-IMEX and C-ETD1 scheme allow for two orders of magnitude larger timesteps than the IMEX and ETD1 scheme, respectively, whereas the C-ETD2RK scheme allows for two times larger timesteps than the ETD2RK scheme. Note, that the (C)-IMEX and (C)-ETD1 integration schemes show similar computational costs for a time step update. Due to an additional intermediate step, the time step updates with the (C)-ETD2RK integration schemes are twice as expensive as with the (C)-IMEX or (C)-ETD1 approaches.

Extension to an explicit temperature coupling
So far, a basic formulation of the PFC model has been considered. However, many extensions have been proposed involving coupling of equation (2) with additional equations. For instance, an extension of the classical PFC model (1) including heat transfer through a temperature field has been proposed in [6,7]. Therein, the (dimensionless) Helmholtz free energy functional x t = 3 · 10 3 t = 2 · 10 3 t = 10 3 (a)  (2) and (6) are shown (b) temperature fields T during crystal growth at t = 10 3 , t = 3 · 10 3 and t = 10 4 for equation (6) are shown (c) time evolution of the free energies F and F β . reads where T (r, t) is the dimensionless temperature with T = 1 at the melting point, and ϑ, γ > 0 are additional parameters. The dynamics are then given by the coupled evolution of T and ψ, reading with M > 0 a parameter corresponding to the thermal diffusivity assumed to be constant. We consider additional initial and boundary conditions for T . In [6], similar concepts to those illustrated in the previous sections have been exploited and are here investigated further. Figure 2 (a)-(c) shows the density and temperature fields of a growing crystal at intermediate time steps with the corresponding free energy evolution. Similar to equation (4) we introduce C 1 , C 2 and report a convergence study of the resulting (C-)IMEX and (C-)ETD1 approaches, in which we compare the residuals of ψ and T , R ψ and R T , respectively, cf. Figure 3. The residuals are evaluated as the discrete L 2 distance from a numerical reference solution, corresponding to a C-IMEX scheme with ∆t = 0.01. As expected, all schemes converge linearly for a decreasing time step size ∆t. Moreover, for a fixed ∆t, the C-IMEX approach gives two orders of magnitude smaller error than the IMEX approach for ψ and almost three orders of magnitude smaller error for T . By fixing an accuracy instead, the C-IMEX approach allows for two orders of magnitude larger timesteps than the IMEX approach, resulting in a further improvement compared to the performances achieved for the numerical integrations as reported in [6]. For the C-ETD1 approach only small performance improvements compared to the ETD1 approach are achieved.

Conclusion
We studied the numerical integration of the governing equations in the PFC model with different time discretizations and featuring a numerical time-stabilization [5]. Also, we showcased by numerical examples that this approach works for both first-order and second-order exponential time integration methods. We outlined the main criteria to determine numerical parameters and investigate the performances for a prototypical simulation, which may be exploited as a benchmark. The applicability to extended PFC models featuring the coupling with additional variables is shown for a model coupling Eq. (2) with heat flow, as reported in [6]. A C-IMEX and a C-ETD1 integration scheme are considered. Especially for the C-IMEX integration scheme, a further improvement is here achieved with respect to a similar approach exploited in [6] thanks to an extended parametrization of the splitting (4).