# Exponential Integrators for Phase-Field Equations using Pseudo-spectral
Methods: A Python Implementation

Elvis do A. Soares<EMAIL_ADDRESS>Amaro G. Barreto Jr  Engenharia de
Processos Químicos e Bioquímicos (EPQB), Escola de Química, Universidade
Federal do Rio de Janeiro, 21941-909, Rio de Janeiro, RJ, Brazil  Frederico W.
Tavares<EMAIL_ADDRESS>Engenharia de Processos Químicos e Bioquímicos
(EPQB), Escola de Química, Universidade Federal do Rio de Janeiro, 21941-909,
Rio de Janeiro, RJ, Brazil  Programa de Engenharia Química, COPPE,
Universidade Federal do Rio de Janeiro, 21941-909, Rio de Janeiro, RJ, Brazil

###### Abstract

In this paper, we implement exponential integrators, specifically Integrating
Factor (IF) and Exponential Time Differencing (ETD) methods, using pseudo-
spectral techniques to solve phase-field equations within a Python framework.
These exponential integrators have showcased robust performance and accuracy
when addressing stiff nonlinear partial differential equations. We compare
these integrators to the well-known implicit-explicit (IMEX) Euler integrators
used in phase-field modeling. The synergy between pseudo-spectral techniques
and exponential integrators yields significant benefits for modeling intricate
systems governed by phase-field dynamics, such as solidification processes and
pattern formation. Our comprehensive Python implementation illustrates the
effectiveness of this combined approach in solving phase-field model
equations. The results obtained from this implementation highlight the
accuracy and computational advantages of the ETD method compared to other
numerical techniques.

##  I Introduction

Phase-field equations [1, 2, 3] are popular approaches for understanding and
modeling complex dynamics, such as phase separation and solidification
processes. The Cahn-Hilliard (CH) model [4, 5, 6] and the Phase-Field Crystal
(PFC) model [7, 8, 9] are the two most well-known phase-field models, with the
former describing phase separation and the latter representing
crystallization. The accurate and efficient numerical solution of these
equations remains an active area of research, with several established methods
tackling the challenges posed by stiff nonlinear partial differential
equations (PDEs) [10]. One of the challenges of phase-field modeling is the
presence of multiple length scales, which can lead to stiff problems that are
computationally expensive to solve. To overcome numerical difficulties,
pseudo-spectral methods have gained attraction because highly accurate and
efficient numerical solutions are obtained.

Pseudo-spectral methods [11, 12] represent a powerful category of numerical
methods capable of solving partial differential equations with periodic
boundary conditions. These techniques deviate by defining representation basis
functions at discrete grid points. This discrete formulation facilitates the
efficient evaluation of specific operators, resulting in substantial
computational speed-up when employing fast algorithms such as the Fast Fourier
Transform (FFT). By discretizing the spatial domain using a Fourier basis and
expressing the solution as a sum of Fourier modes, pseudo-spectral methods
achieve high accuracy while minimizing computational costs. These methods have
been successfully applied to a diverse range of problems in physics and
chemistry, including Burgers, Navier-Stokes [13], Kardar-Parisi-Zhang [14],
Allen-Cahn [15], Cahn-Hilliard [16, 17], and Phase-Field-Crystal equations
[18, 19]. As pseudo-spectral methods have exhibited strong performance in this
variety of applications, it is crucial to explore how their combination with
advanced time integration schemes can further enhance the solution of stiff
nonlinear PDEs.

Recently, three robust numerical methods have gained prominence for solving
stiff nonlinear PDEs: Implicit-Explicit (IMEX), Integrating Factor (IF) [20,
21, 22, 23], and Exponential Time Differencing (ETD) methods [24, 25, 26].
These exponential integrators have showcased robust performance and accuracy
when addressing stiff nonlinear partial differential equations. These methods
present an appealing alternative to the widely-used IMEX Euler integrators in
phase-field modeling. By capitalizing on the synergy between pseudo-spectral
techniques and exponential integrators, we aim to emphasize the significant
benefits and improvements in modeling complex systems governed by phase-field
dynamics.

Building upon the demonstrated advantages of combining pseudo-spectral
techniques with exponential integrators, this paper probes a comprehensive
Python implementation of these integrators, showcasing their effectiveness in
solving phase-field model equations and emphasizing the potential of Python
usage for advancing phase-field modeling. The results obtained from our Python
implementation stress the accuracy and computational efficiency of the IMEX
and ETD methods, further highlighting the potential of this approach in
advancing phase-field modeling research. The codes presented in this work are
available online (https://github.com/elvissoares/spectralETD) under a GPL
license.

The paper is organized as follows: In Section II, a brief review of the
theoretical formalism for phase-field models, pseudo-spectral methods, and
time integration schemes is presented. The numerical implementation is
detailed in Section III. Results of the numerical calculations for phase field
models are provided in Section IV. Section V concludes this paper.

##  II Formalism

###  II.1 Phase-Field Models with Conserved Order Parameters

Phase-Field models [27] are based on a mathematical concept called a phase
field, $\eta({\bf\it r},t)$, an order parameter that is a continuous function
that describes the spatial distribution of different phases or microstructures
in a material.

Here we concentrate on the case where $\eta({\bf\it r},t)$ is a conserved
order parameter. The conservation equation for $\eta({\bf\it r},t)$ is
described by $\partialderivative*{\eta}{t}+\divergence{{\bf\it j}_{\eta}}=0$,
with the flux density given by ${\bf\it j}_{\eta}=-M\nabla(\delta
F[\eta]/\delta\eta)$, where $F[\eta]$ is the free-energy functional of the
phase field. The dynamics of the system can be described by

$\displaystyle\partialderivative{\eta}{t}=\divergence\left[M(\eta)\nabla\left(\functionalderivative{F[\eta]}{\eta}\right)\right],$
(1)

where $M(\eta)$ is a positive-defined mobility coefficient that can depend on
$\eta$.

The Cahn-Hilliard (CH) equation is a phase-field model that describes the
dynamics of phase separation and coarsening process in two-component mixtures,
and it has been applied to various other fields, such as material science,
biology, and fluid dynamics. It was introduced in 1958 by John W. Cahn and
John E. Hilliard [4] as a mathematical model to understand the thermodynamics
and kinetics of phase separation in alloys at a mesoscopic scale, between the
atomic and macroscopic scales. The CH free-energy functional is defined as

$\displaystyle F[\eta({\bf\it r})]=\int_{V}\differential{{\bf\it
r}}\left[\frac{\kappa}{2}\left(\nabla\eta\right)^{2}+f(\eta)\right],$ (2)

where $\eta({\bf\it r})$ is the concentration field of one phase, $\kappa$ is
a positive constant that controls the interfacial energy between the two
phases, and the term $f(\eta)=W\eta^{2}(1-\eta)^{2}$ defines the double-well
free-energy density function, which describes the energetics of the phase
separation process. The Eq. (1) leads to the dynamical equation given by

$\displaystyle\partialderivative{\eta}{t}=M\laplacian\left[-\kappa\laplacian\eta+2W(\eta-3\eta^{2}+2\eta^{3})\right],$
(3)

with the $\kappa\laplacian\eta$ term accounting for the interfacial energy.
This term tends to smooth out the concentration field and promote the
formation of well-defined interfaces between the two phases. Therefore, the CH
equation models the time evolution of the order parameter $\eta$ as a result
of the competition between the driving force for phase separation, which comes
from the free-energy density function $f(\eta)$, and the energy penalty
$\kappa\laplacian\eta$ associated with the formation of interfaces between the
phases.

The phase field crystal (PFC) equation describes crystallization growth, where
atomic- and microscales are coupled. Like conventional phase field models, the
PFC theory incorporates a free energy functional and averages over rapid
temporal fluctuations, thereby yielding a time scale of evolution that is on
the order of diffusion rather than atomic vibrations. Besides that, the PFC
model does not average over atomic distances, leading to the formation of
equilibrium patterns that are representative of the crystalline structures.

Here we focus on the model from Elder and Grant presented in Ref. [7]. In this
PFC model, the free-energy functional is written as

$\displaystyle F[\eta({\bf\it r})]=\int_{V}\differential{{\bf\it
r}}\left[\frac{1}{2}\eta\left(1+\laplacian\right)^{2}\eta+\frac{1}{4}\eta^{2}(2r+\eta^{2})\right],$
(4)

with $r$ being a constant proportional to the temperature deviation from the
melting point, and the term $\frac{1}{2}\eta\left(1+\laplacian\right)^{2}\eta$
reproduces the first-order peak of the structure factor for the liquid phase.

Using the functional derivative of this free-energy in Eq. (1), we get

$\displaystyle\partialderivative{\eta}{t}=M\laplacian\quantity[\nabla^{4}\eta+2\laplacian\eta+(1+r)\eta+\eta^{3}],$
(5)

where the first two terms account for the formation of the interface between
the liquid-solid introducing the periodic pattern in the solid phase,
respectively. The last two terms in Eq. (5) constitute a double-well free
energy between the liquid and solid phases.

All the phase field models with conserved order parameters have conservative
dynamics that minimize the total free-energy functional, $F[\eta({\bf\it
r},t)]$, of the system. In fact, the dynamics described by Eq. (1) leads to
the free-energy evolution given by

$\displaystyle\derivative{F}{t}=-\int_{V}\differential{{\bf\it
r}}M(\eta)\quantity[\gradient(\functionalderivative{F[\eta]}{\eta})]^{2}\geq
0,$ (6)

such that the free energy can just decrease or remain constant during the time
evolution of $\eta({\bf\it r},t)$.

###  II.2 Pseudo-Spectral Version of Conservative Equations

The dynamics of any field $\eta({\bf\it r},t)$ described by the Eq. (1) can be
written as the sum of a set of linear terms and another set of non-linear
terms in the form

$\displaystyle\partialderivative{\eta}{t}=\sum_{\alpha=0}^{A}\mathcal{L}_{\alpha}\nabla^{\alpha}\eta+\sum_{\beta=0}^{B}\mathcal{N}_{\beta}\nabla^{\beta}f_{\beta}(\eta),$
(7)

where $\eta$ is the field, $\nabla$ is the spatial coordinate derivative
operator, $t$ is time, $\mathcal{L}_{\alpha}$ and $\mathcal{N}_{\beta}$ are
numerical coefficients of the linear and non-linear terms, respectively, and
$f_{\beta}(\eta)$ are some non-linear functions of the field $\eta({\bf\it
r},t)$. Here, $\alpha$ and $\beta$ are nonnegative integer indexes for the
$A+1$ linear terms and $B+1$ non-linear terms, respectively. On conservative
systems, we expect that $\alpha\geq 2$.

We can decompose the field $\eta$ into independent normal modes by setting the
direct Fourier transform as

$\displaystyle\widetilde{\eta}_{{\bf\it k}}(t)\equiv\mathcal{F}\\{\eta({\bf\it
r},t)\\}=\int_{L^{d}}\differential{{\bf\it r}}\eta({\bf\it r},t)e^{-i{\bf\it
k}\cdot{\bf\it r}},$ (8)

and the inverse Fourier transform as

$\displaystyle\eta({\bf\it
r},t)\equiv\mathcal{F}^{-1}\\{\widetilde{\eta}_{{\bf\it
k}}(t)\\}=\frac{1}{N^{d}}\sum_{{\bf\it k}}\widetilde{\eta}_{{\bf\it
k}}(t)e^{i{\bf\it k}\cdot{\bf\it r}},$ (9)

with ${\bf\it k}=\sum_{i=1}^{d}k_{i}\hat{{\bf\it u}}_{i}$ being the wavenumber
vector in the Fourier space in $d$ dimensions, and $\mathcal{F}$ and
$\mathcal{F}^{-1}$ represent the direct and inverse Fourier transforms,
respectively. The complete Fourier space can be defined in each direction by
the wavenumbers given by $k_{i}=\\{0,2\pi/L,4\pi/L,\ldots,2\pi N/L\\}$, where
$L$ and $N$ are the length of the box and the number of gridpoints in the
$i$-direction, respectively.

Any spatial derivative of the field $\eta({\bf\it r},t)$ can be calculated by
the Fourier transform as $\nabla^{\alpha}\eta=\sum_{{\bf\it k}}(i{\bf\it
k})^{\alpha}\widetilde{\eta}_{{\bf\it k}}e^{i{\bf\it k}\cdot{\bf\it r}}$.
Similarly, we have $\nabla^{\beta}f_{\beta}(\eta)=\sum_{{\bf\it k}}(i{\bf\it
k})^{\beta}\mathcal{F}\quantity{f_{\beta}(\eta)}_{{\bf\it k}}e^{i{\bf\it
k}\cdot{\bf\it r}}$. Collecting the Fourier coefficients, we have

$\displaystyle\partialderivative{\widetilde{\eta}_{{\bf\it
k}}}{t}=\sum_{\alpha=0}\mathcal{L}_{\alpha}(i{\bf\it
k})^{\alpha}\widetilde{\eta}_{{\bf\it
k}}+\sum_{\beta=0}\mathcal{N}_{\beta}(i{\bf\it
k})^{\beta}\mathcal{F}\quantity{f_{\beta}(\eta)}_{{\bf\it k}},$ (10)

and defining the auxiliary linear and non-linear marching operators,
respectively, in the form $\widetilde{\mathcal{L}}_{{\bf\it
k}}=\sum_{\alpha=0}\mathcal{L}_{\alpha}(i{\bf\it k})^{\alpha},$ and
$\widetilde{\mathcal{N}}_{{\bf\it
k}}(\eta(t))=\sum_{\beta=0}\mathcal{N}_{\beta}(i{\bf\it
k})^{\beta}\mathcal{F}\quantity{f_{\beta}(\eta)}_{{\bf\it k}}$, we can
simplify the FT dynamical equation as

$\displaystyle\partialderivative{\widetilde{\eta}_{{\bf\it
k}}}{t}=\widetilde{\mathcal{L}}_{{\bf\it k}}\widetilde{\eta}_{{\bf\it
k}}+\widetilde{\mathcal{N}}_{{\bf\it k}}(\eta(t)).$ (11)

Note that the linear marching operator is time-independent while the non-
linear marching operator is time-dependent because the non-linear functions
$f_{\beta}$ are dependent on $\eta(t)$.

If there aren’t any non-linear terms, _i.e._ , $f_{\beta}=0$, the equation has
an analytical solution in the form $\widetilde{\eta}_{{\bf\it
k}}(t)=\widetilde{\eta}_{{\bf\it k}}(0)e^{t\widetilde{\mathcal{L}}_{{\bf\it
k}}}$. Therefore, the product $\widetilde{\eta}_{{\bf\it
k}}(t)e^{-t\widetilde{\mathcal{L}}_{{\bf\it k}}}$ is a time-invariant of the
pure linear model. In this manner, the exponential term
$e^{-t\widetilde{\mathcal{L}}_{{\bf\it k}}}$ can be used as an integrand
factor to solve the more general problem, Eq. (11). Multiplying the Eq. (11)
by the integrand factor $e^{-t\widetilde{\mathcal{L}}_{{\bf\it k}}}$, we can
identify

$\displaystyle\partialderivative{t}\quantity(\widetilde{\eta}_{{\bf\it
k}}(t)e^{-t\widetilde{\mathcal{L}}_{{\bf\it
k}}})=\widetilde{\mathcal{N}}_{{\bf\it
k}}(\eta(t))e^{-t\widetilde{\mathcal{L}}_{{\bf\it k}}},$ (12)

without any loss of generality. Integrating Eq. (12) from $t_{0}$ to $t$, we
get

$\displaystyle\widetilde{\eta}_{{\bf\it
k}}(t)=\quantity[\widetilde{\eta}_{{\bf\it
k}}(t_{0})+\int_{0}^{t-t_{0}}\widetilde{\mathcal{N}}_{{\bf\it
k}}(\eta(\tau+t_{0}))e^{-\tau\widetilde{\mathcal{L}}_{{\bf\it
k}}}\differential{\tau}]e^{(t-t_{0})\widetilde{\mathcal{L}}_{{\bf\it k}}}.$
(13)

which is an analytical relation without any approximations.

###  II.3 Time Integration Methods

####  II.3.1 Implicit-Explicit Euler Method: IMEX

The IMEX (IM = implicit, EX = explicit) method [17] can be used to approximate
the time integration of Eq. (11). The method can be obtained by discretizing
Eq. (11) from $t_{n}$ until $t_{n+1}=t_{n}+h$, using an implicit step for the
linear term and an explicit step for the nonlinear term, as follows

$\displaystyle\frac{\widetilde{\eta}_{{\bf\it
k}}(t_{n+1})-\widetilde{\eta}_{{\bf\it
k}}(t_{n})}{h}=\widetilde{\mathcal{L}}_{{\bf\it k}}\widetilde{\eta}_{{\bf\it
k}}(t_{n+1})+\widetilde{\mathcal{N}}_{{\bf\it k}}(\eta(t_{n})),$ (14)

where $t_{n+1}=t_{n}+h$. The Eq. (14) can be solved for
$\widetilde{\eta}_{{\bf\it k}}(t_{n+1})$ such that

$\displaystyle\widetilde{\eta}_{{\bf\it
k}}(t_{n+1})=\frac{\widetilde{\eta}_{{\bf\it
k}}(t_{n})+h\widetilde{\mathcal{N}}_{{\bf\it
k}}(\eta(t_{n}))}{1-h\widetilde{\mathcal{L}}_{{\bf\it k}}}.$ (15)

The implicit scheme is stable, so larger time steps may be used. But there is
a necessary condition of $h|\widetilde{\mathcal{L}}_{{\bf\it k}}|<1.0$ to make
the convergence.

####  II.3.2 Integrating Factor Method: IF

The IF method [20, 21, 22, 23] can be obtained by discretizing the Eq. (12) in
time from $t_{n}$ until $t_{n+1}=t_{n}+h$, using an explicit step for the
nonlinear source term in the following form

$\displaystyle\frac{\widetilde{\eta}_{{\bf\it
k}}(t_{n+1})e^{-t_{n+1}\widetilde{\mathcal{L}}_{{\bf\it
k}}}-\widetilde{\eta}_{{\bf\it
k}}(t_{n})e^{-t_{n}\widetilde{\mathcal{L}}_{{\bf\it
k}}}}{h}=\widetilde{\mathcal{N}}_{{\bf\it
k}}(\eta(t_{n}))e^{-t_{n}\widetilde{\mathcal{L}}_{{\bf\it k}}},$ (16)

such that

$\displaystyle\widetilde{\eta}_{{\bf\it k}}(t_{n+1})=\widetilde{\eta}_{{\bf\it
k}}(t_{n})e^{h\widetilde{\mathcal{L}}_{{\bf\it
k}}}+h\widetilde{\mathcal{N}}_{{\bf\it
k}}(\eta(t_{n}))e^{h\widetilde{\mathcal{L}}_{{\bf\it k}}}.$ (17)

If we multiply Eq. (17) by $e^{-h\widetilde{\mathcal{L}}_{{\bf\it k}}}$, and
take the limit $h|\widetilde{\mathcal{L}}_{{\bf\it k}}|\ll 1$ such that
$e^{-h\widetilde{\mathcal{L}}_{{\bf\it k}}}\approx
1-h\widetilde{\mathcal{L}}_{{\bf\it k}}$, we obtain the IMEX method.

####  II.3.3 Exponential Time Differencing Method: ETD

The ETD method [28, 25] consists of discretizing Eq. (13) in time from $t_{n}$
until $t_{n+1}=t_{n}+h$. Using a simple constant approximation of the non-
linear operator in the integral $\int_{0}^{h}\widetilde{\mathcal{N}}_{{\bf\it
k}}(\eta(\tau+t_{n}))e^{-\tau\widetilde{\mathcal{L}}_{{\bf\it
k}}}\differential{\tau}\approx\widetilde{\mathcal{N}}_{{\bf\it
k}}(\eta(t_{n}))[1-e^{-h\widetilde{\mathcal{L}}_{{\bf\it
k}}}]/\widetilde{\mathcal{L}}_{{\bf\it k}}$, such that

$\displaystyle\widetilde{\eta}_{{\bf\it k}}(t_{n+1})=\widetilde{\eta}_{{\bf\it
k}}(t_{n})e^{h\widetilde{\mathcal{L}}_{{\bf\it
k}}}+h\widetilde{\mathcal{N}}_{{\bf\it
k}}(\eta(t_{n}))\quantity[\frac{e^{h\widetilde{\mathcal{L}}_{{\bf\it
k}}}-1}{h\widetilde{\mathcal{L}}_{{\bf\it k}}}].$ (18)

Note that, multiplying Eq. (18) by $e^{-h\widetilde{\mathcal{L}}_{{\bf\it
k}}}$ and taking $|\widetilde{\mathcal{L}}_{{\bf\it k}}|\ll 1$, we obtain the
IMEX method.

The function $(e^{x}-1)/x$, despite being continuous, presents undesirable
behavior near $x=0$ in a computational context due to cancellation errors.
This can be bypassed by using the appropriate limit of the function at $x=0$.

##  III Numerical Implementation in Python

We have implemented all the time integration schemes combined with the pseudo-
spectral method, as described in the previous section, in high-level Python
code.

For the 1D geometry, we comment on the cases of the advection-diffusion
equation and the Burgers’ equation in Appendix A.

In the examples considered here, we consider a 2D geometry as our physical
space that consists of a structured uniform square $[0,16\pi]^{2}$, where each
direction is divided into $N=2^{M}$ uniform intervals, where $M$ is a positive
integer. The computational mesh is represented in Python as

⬇

1import numpy as np

2# Size of the system

3N = 2**8 # 2**8 = 256

4L = 16*np.pi

5x = np.linspace(0,L,N)

6dx = x[1]-x[0]

The system was evolved with a stepsize of $h=0.01$, and simulation was carried
out up to 150 000 time steps. At each 100th frame, we output and store the
result in an array.

⬇

1# The time step definition

2h = 0.01

3T = 1500

4Nsteps = int(T/h)

5dframes = 1.0 # time step to output

6Nframes = int(T/dframes) #frames to the output

7nframes = Nsteps//Nframes

8# The array of outputs

9n = np.empty((Nframes,N,N), dtype=np.float32)

Our pseudo-spectral Python solver extensively uses the NumPy [29] package. The
FFT module comes from Scipy [30] package. The fourier space can be defined by

⬇

1# The Fourier variables

2from scipy.fft import fft2, ifft2

3n_k = np.empty((N,N), dtype=np.complex64)

4kx = np.fft.fftfreq(N, d=dx)*2*np.pi

5k = np.array(np.meshgrid(kx , kx ,indexing =’ij’), dtype=np.float32)

6k2 = np.sum(k*k,axis=0, dtype=np.float32)

To avoid any aliasing problem with non-linear terms due to the discrete nature
of numerical simulations, we have used the _2/3 rule_ dealising technique
[31]. The aliasing happens when the signal or function is not adequately
sampled, violating the Nyquist-Shannon sampling theorem. Here, we introduce
the dealising matrix as $\xi^{\text{dealias}}_{{\bf\it
k}}=\Pi_{i}^{d}\Theta\quantity(k_{i}^{\text{cut}}-k_{i})$, where
$k_{i}^{\text{cut}}=(2/3)k_{i}^{\text{max}}$ is the Nyquist cutoff frequency
in the $i$th direction and $\Theta(x)$ is the Heaviside step function.

⬇

1kmax_dealias = kx.max()*2.0/3.0 # The Nyquist mode

2# Dealising matrix

3dealias = np.array((np.abs(k[0]) < kmax_dealias )*(np.abs(k[1]) <
kmax_dealias ),dtype =bool)

The linear and nonlinear marching operators can be defined by

⬇

1# The linear terms of PDE

2Loperator_k = (...) # some function of k only

3# The non-linear terms of PDE

4def Noperator_func(n):

5 return (...) # some function of n and k

After that, the time integration schemes can be chosen using the following
code

⬇

1# Defining the time marching operators arrays

2# can be calculated once

3if method == ’IMEX’:

4 Tlinear_k = 1.0/(1.0-h*Loperator_k)

5 Tnon_k = dealias*h/(1.0-h*Loperator_k)

6elif method == ’IF’:

7 Tlinear_k = np.exp(h*Loperator_k)

8 Tnon_k = dealias*h*Tlinear_k

9elif method == ’ETD’:

10 Tlinear_k = np.exp(h*Loperator_k)

11 def myexp(x):

12 if x == 1: return 1.0

13 else: return (x-1.0)/np.log(x)

14 vmyexp = np.vectorize(myexp) # vectorize myexp

15 Tnon_k = dealias*h*vmyexp(Tlinear_k)

16else: print(’ERROR: Undefined Integrator’)

Note that the operator can be calculated once and before the time integration
loop.

Finally, the loop of time integration is

⬇

1n[0] = (...) # some initial condition

2Noperator_k = n_k.copy() # auxiliary array

3nn = n[0].copy() # auxiliary array

4n_k[:] = fft2(n[0]) # FT initial condition

5# time evolution loop

6for i in range(1,Nsteps):

7 # calculate the nonlinear operator

8 Noperator_k[:] = Noperator_func(nn)

9 # updating in time

10 n_k[:] = n_k*Tlinear_k + Noperator_k*Tnon_k

11 # IFT to next step

12 nn[:] = ifft2(n_k).real

13 # test to output

14 if (i % nframes) == 0: n[i//nframes] = nn

##  IV Results and Discussion

###  IV.1 2D Cahn-Hilliard Equation

Here we consider the CH equation with the parameters set as $W=1$,
$\kappa=0.1$, and $M=1$. The initial profile is a homogeneous field with value
$\eta_{0}$ and a Gaussian perturbation of intensity 0.02 at each gridpoint.

⬇

1# Cahn-Hilliard model constants

2W = 1.0

3M = 1.0 # mobility

4kappa = 0.1 #gradient coeficient

5# Initial condition

6rng = np.random.default_rng(12345) # the seed

7noise = 0.02

8n0 = 0.5

9n[0] = n0 +noise*rng.standard_normal(n[0].shape)

From the pseudo-spectral method, the CH equation of motion, Eq. (3), can be
re-written as Eq. (11) with the time-marching operators given by

$\displaystyle\widetilde{\mathcal{L}}_{{\bf\it k}}$ $\displaystyle=-M(\kappa
k^{4}+2Wk^{2}),$ (19) $\displaystyle\widetilde{\mathcal{N}}_{{\bf\it
k}}(\eta)$
$\displaystyle=-2MWk^{2}\mathcal{F}\quantity{-3\eta^{2}+2\eta^{3}}_{{\bf\it
k}},$ (20)

with the following Python code:

⬇

1# The linear terms of PDE

2Loperator_k = -M*(kappa*k2**2+2*W*k2)

3# The non-linear terms of PDE

4def Noperator_func(n):

5 return -2*M*W*k2*fft2(-3*n**2+2*n**3)

The whole Python code to solve the CH equation with exponential integrators
and pseudo-spectral methods is presented in Appendix B.

Figure 1 shows the result of the CH equation using the pseudo-spectral method
with the ETD integration scheme for three different initial condition
$\eta_{0}=0.4$, 0.5 and 0.6 (see Movies S1, S2, and S3 for more details). The
first row presents the dynamics of the nucleation and growth of bubbles of the
red phase inside the continuous blue phase. The second row represents the
dynamics of spinodal decomposition of the two phases. And the third row
presents the nucleation dynamics of the blue phase surrounded by the
continuous red phase.

Figure 1: The time evolution of the field profile $\eta({\bf\it r},t)$
calculated from the CH equation. The different rows represent different
initial conditions: $\eta_{0}=0.4$ (see Movie S1), $\eta_{0}=0.5$ (see Movie
S2), and $\eta_{0}=0.6$ (see Movie S3). The columns represents the snapshots
in different instant of time $t$. The system was evolved with the ETD scheme
and a stepsize of $h=0.01$. The initial profile is a homogeneous field with
$\eta_{0}$ and a Gaussian perturbation of intensity 0.02, _i.e._ ,
$\eta(t=0)=\eta_{0}+0.02\mathcal{N}(0,1)$.

Figure 2 illustrates the temporal evolution of the CH free energy, as defined
by Eq. (2), during spinodal decomposition ($\eta_{0}=0.5$) for the three
distinct time integration schemes. It can be observed intense decay in the
free-energy during the initial times up to $t=100$. Following that, the system
relaxes at a slower rate for approximately 1400 units of time. The various
integration schemes display minor variations in the free energy decay curve
over the course of the process, although they converge to a similar energy
value towards the end of the numerical calculation.

Figure 2: Free energy as a function of time to the CH equation with the initial condition as a homogeneous field with $\eta_{0}=0.5$ and a Gaussian perturbation of intensity 0.02 with integrating using different methods with a timestep of $h=0.01$. The solid line represents the solution with the ETD method, the dashed line with the IF method, and the dotted line with the IMEX method. Table 1: The $L_{2}$-error of the numerical solution of the CH equation for the spinodal decomposition condition with $\eta_{0}=0.5$ at $t=100$. $h$ | IMEX | IF | ETD
---|---|---|---
0.0002 | $6.4\times 10^{-4}$ | $9.4\times 10^{-4}$ | $6.4\times 10^{-4}$
0.001 | $7.4\times 10^{-4}$ | $1.3\times 10^{-3}$ | $6.1\times 10^{-4}$
0.002 | $7.5\times 10^{-4}$ | $1.7\times 10^{-3}$ | $7.4\times 10^{-4}$
0.01 | $1.4\times 10^{-3}$ | $2.2\times 10^{-3}$ | $1.3\times 10^{-3}$
0.02 | $1.8\times 10^{-3}$ | $2.4\times 10^{-3}$ | $1.4\times 10^{-3}$
0.1 | $2.0\times 10^{-3}$ | $2.4\times 10^{-3}$ | $2.0\times 10^{-3}$
0.2 | $2.2\times 10^{-3}$ | $2.2\times 10^{-3}$ | $2.2\times 10^{-3}$
1.0 | $2.3\times 10^{-3}$ | $1.9\times 10^{-3}$ | $2.4\times 10^{-3}$

The $L_{2}$-error is defined as the absolute error given by
$\norm{\eta-\eta_{r}}_{2}/N^{2}$ where $\eta$ is the numerical solution with
determined stepsize $h$ and determined scheme, and $\eta_{r}$ is the numerical
solution with the reference stepsize of $10^{-4}$ and ETD scheme. The
normalization $N^{2}$ is the total number of gridpoints. Table 1 presents the
$L_{2}$-error of the three different time integration schemes. The ETD method
comes up as the most accurate when compared to the IF method. However, the
simple IMEX method remains competitive even when taking larger time steps into
account.

###  IV.2 2D Phase Field Crystal Equation

We consider the PFC equation with the parameters set as $r=-0.25$ and $M=1.0$.
The initial profile is a homogeneous field $\eta_{0}$ with a Gaussian
perturbation of intensity $0.02|\eta_{0}|$ at each gridpoint.

⬇

1# PFC model constants

2r = -0.25

3M = 1.0 # mobility

4# Initial condition

5rng = np.random.default_rng(12345)

6n0 = -0.285

7noise = 0.02*np.abs(n0)

8n[0] = n0 +noise*rng.standard_normal(n[0].shape)

From the pseudo-spectral method, the PFC equation can be re-written as Eq.
(11) with the time-marching operators given by

$\displaystyle\widetilde{\mathcal{L}}_{{\bf\it k}}$
$\displaystyle=-Mk^{2}(k^{4}-2k^{2}+1+r),$ (21)
$\displaystyle\widetilde{\mathcal{N}}_{{\bf\it k}}(\eta)$
$\displaystyle=-Mk^{2}\mathcal{F}\quantity{\eta^{3}}_{{\bf\it k}},$ (22)

with the following Python code:

⬇

1# The linear terms of PDE

2Loperator_k = -M*k2*(k2**2-2*k2+1+r)

3# The non-linear terms of PDE (with dealising)

4def Noperator_func(n):

5 return -(k2*M*fft2(n**3))

Again, the whole Python code to solve the PFC equation with exponential
integrators and pseudo-spectral methods is presented in Appendix B.

Figure 3 presents the time evolution of the field profile calculated from the
PFC equation using the pseudo-spectral method with the ETD integration scheme.
The first line presents the snapshots of the temporal evolution of a system
with an initial condition of $\eta_{0}=-0.085$, which leads to the formation
of a lamellar phase. The second line presents the snapshots of the system’s
temporal evolution with $\eta_{0}=-0.285$, which results in an ordered
crystalline phase. (see Movies S4 and S5 for more details).

Figure 3: The time evolution of the field profile $\eta({\bf\it r},t)$
calculated from the PFC equation. The different rows represent the two
different initial conditions characterized by $\eta_{0}=-0.085$ (see Movie S4)
and $\eta_{0}=-0.085$ (see Movie S5). The columns represent the snapshots of
different time instant $t$. The system was evolved with an ETD scheme and
stepsize of $h=0.01$. The initial profile is a homogeneous field with
$\eta_{0}$ and a Gaussian perturbation of intensity $0.02|\eta_{0}|$, _i.e._ ,
$\eta(t=0)=\eta_{0}+0.02|\eta_{0}|\mathcal{N}(0,1)$.

Figure 4 represents the temporal evolution of the free energy for the system
generating a lamellar phase derived from the PFC equation. A significant
topological transition occurs around $t=10$, leading to the initial formation
of a lamellar phase. Between the time interval from $t=100$ to $t=300$, the
lamellar phase forms throughout the entire system. The free energy remains
constant after this period until the final time step.

Figure 4: Free energy as a function of time from the PFC equation with the
initial condition as a homogeneous field with $\eta_{0}=-0.085$ and a Gaussian
perturbation of intensity $0.02|\eta_{0}|$.

Figure 5 presents the free-energy as a function of time during the crystal
formation on the PFC dynamics. There is a sharp decay of the free-energy until
$t=10$ due to the decay of the high energy modes. Between $t=10$ and 700, the
system experiences a slight variation in free energy, with minimal movements
in the phase field configuration. However, near $t=750$, there is a rapid
decrease in free energy during the topological transition to form the ordered
crystalline lattice. Indeed, $t=750$ is the time required for the system to
form the crystal and is associated with the system’s crystalline nucleation
rate. It is interesting to note that this jump in free energy only occurs
during the formation of crystalline phases and does not occur during the
formation of the lamellar phase, according to Fig. 4.

Figure 5: Free energy as a function of time from the PFC equation with the initial condition as a homogeneous field with $\eta_{0}=-0.285$ and a Gaussian perturbation of intensity $0.02|\eta_{0}|$.  Table 2: The $L_{2}$-error of the numerical solution of PFC equation for the crystallization condition with $\eta_{0}=-0.285$ at $t=750$. $h$ | IMEX | IF | ETD
---|---|---|---
0.0002 | $1.5\times 10^{-6}$ | $1.5\times 10^{-6}$ | $1.8\times 10^{-6}$
0.001 | $5.5\times 10^{-6}$ | $5.4\times 10^{-6}$ | $7.0\times 10^{-6}$
0.002 | $5.5\times 10^{-6}$ | $2.3\times 10^{-6}$ | $5.5\times 10^{-6}$
0.01 | $5.1\times 10^{-6}$ | $1.2\times 10^{-5}$ | $4.8\times 10^{-6}$
0.02 | $5.3\times 10^{-6}$ | $3.1\times 10^{-5}$ | $4.8\times 10^{-6}$
0.1 | $9.5\times 10^{-6}$ | $1.8\times 10^{-4}$ | $6.9\times 10^{-6}$
0.2 | $1.4\times 10^{-5}$ | $2.3\times 10^{-4}$ | $9.4\times 10^{-6}$
1.0 | $3.0\times 10^{-5}$ | $2.3\times 10^{-4}$ | $2.3\times 10^{-5}$

Table 2 presents the $L_{2}$-error of the three different time integration
schemes in the case of the PFC equation for crystalization. The ETD method
comes up as the most accurate when compared to the other two methods. The ETD
method remains very precise even when taking larger time steps into account.
It is worth noting that the IF method leads to a different nucleation rate
when compared to the other two methods. Therefore, avoiding the IF method when
dealing with PFC equations is interesting.

##  V Conclusions

A comprehensive investigation has been conducted, focusing on three distinct
time integration schemes (IMEX, IF, and ETD) in combination with pseudo-
spectral methods, which are designed to lead with stiff or highly oscillatory
nonlinear PDEs, thereby highlighting the potential utility of these
methodologies in addressing complex scientific problems.

The ETD method demonstrates its superiority in accuracy and enhanced stability
compared to the IF method, solidifying its position as a preferred choice for
handling complex numerical problems. Moreover, avoiding the IF method when
addressing PFC equations is important, as it results in a distinct nucleation
rate compared to the other two methods. However, it is worth noting that the
simple IMEX method maintains its competitiveness. This suggests that,
depending on the specific requirements and constraints of a given problem, the
IMEX method may still be a viable alternative for first applications. After
all, the ETD method exhibits enhanced stability compared to other methods,
further contributing to its appeal as a robust and reliable choice for
tackling various phase-field modeling scenarios. In summary, weighing the
trade-offs between accuracy, stability, and computational efficiency is
essential when selecting the most appropriate method for a particular phase-
field modeling scenario.

The combination of pseudo-spectral techniques and exponential integrators
showcases significant benefits for modeling complex systems governed by phase-
field dynamics, including solidification processes and pattern formation. Our
comprehensive Python implementation shows the effectiveness of this combined
approach in solving phase-field model equations, highlighting the accuracy and
computational advantages of the ETD method compared to other numerical
approaches.

In conclusion, this work underscores the potential of leveraging exponential
integrators and pseudo-spectral techniques to advance phase-field modeling
research and provide accurate and efficient solutions for complex scientific
problems. We aim to expand this methodology by improving ETD via Runge-Kutta
schemes, as already discussed in Ref. [25] for other equations. We are also
interested in examining its applicability in the context of stochastic partial
differential equations.

## Acknowledgments

The authors thank Petrobras and Shell, which provided financial support
through the Research, Development, and Innovation Investment Clause, in
collaboration with the Brazilian National Agency of Petroleum, Natural Gas,
and Biofuels (ANP, Brazil). Additionally, this research was partially funded
by CNPq, CAPES, and FAPERJ.

## The data and code availability

The data and code that support the findings of this study are available in the
author’s Github repository: https://github.com/elvissoares/spectralETD.

## Supplementary Material

Supplementary material for this article is available at: link

  * •

_Movie S1_. Time evolution of 2D CH equation ($W=1$, $\kappa=0.1$ and $M=1$)
with $\eta_{0}=0.4$ (red phase nucleation-growth-ripening) using pseudo-
spectral method with ETD scheme.

  * •

_Movie S2_. Time evolution of 2D CH equation ($W=1$, $\kappa=0.1$ and $M=1$)
with $\eta_{0}=0.5$ (spinodal decomposition) using pseudo-spectral method with
ETD scheme.

  * •

_Movie S3_. Time evolution of 2D CH equation ($W=1$, $\kappa=0.1$ and $M=1$)
with $\eta_{0}=0.6$ (blue phase nucleation-growth-ripening) using pseudo-
spectral method with ETD scheme.

  * •

_Movie S4_. Time evolution of 2D PFC equation ($r=-0.25$ and $M=1$) with
$\eta_{0}=-0.085$ (lamellar phase) using pseudo-spectral method with ETD
scheme.

  * •

_Movie S5_. Time evolution of 2D PFC equation ($r=-0.25$ and $M=1$) with
$\eta_{0}=-0.285$ (crystal phase) using pseudo-spectral method with ETD
scheme.

  * •

_Movie S6_. Time evolution of 1D advection-diffusion equation ($u=5$ and
$D=0.01$) using the analytical solution in Fourier space.

  * •

_Movie S7_. Time evolution of 1D Burgers’ equation ($\nu=0.001$) using pseudo-
spectral method with ETD scheme.

##  Appendix A Other Examples

###  A.1 1D Advection-Diffusion Equation

The advection-diffusion equation is a fundamental partial differential
equation that describes the transport of a scalar quantity (e.g., temperature,
concentration) in a fluid flow. It combines the effects of advection,
representing the transport of the scalar quantity by the fluid flow, and
diffusion, which represents the random movement of the scalar quantity due to
molecular motion. The one-dimensional advection-diffusion equation is written
as

$\displaystyle\frac{\partial\eta}{\partial t}=-u\frac{\partial\eta}{\partial
x}+D\frac{\partial^{2}\eta}{\partial x^{2}},$ (23)

with $u$ being the velocity constant, $D$ being the diffusion coefficient, and
$\eta$ represents the concentration field. The time derivative of the Fourier
Transform $\widetilde{\eta}_{{\bf\it k}}$ is

$\displaystyle\partialderivative{t}\widetilde{\eta}_{{\bf\it
k}}(t)=-(iuk+Dk^{2})\widetilde{\eta}_{{\bf\it k}}(t).$ (24)

This case corresponds to the pseudo-spectral equation Eq. (11) with the
following operators given by

$\displaystyle\widetilde{\mathcal{L}}_{{\bf\it
k}}=-(iuk+Dk^{2}),\quad\text{and}\quad\widetilde{\mathcal{N}}_{{\bf\it k}}=0,$
(25)

such that the time evolution has an analytical solution in the form

$\displaystyle\widetilde{\eta}_{k}(t)=\widetilde{\eta}_{k}(0)e^{-(iuk+Dk^{2})t}.$
(26)

As an example, we set $u=5$ and $D=0.01$. The size of the system is $L=2\pi$
with the number of gridpoints $N=2^{12}=4096$. The initial profile is a top-
hat function with a length of $l=0.2$ and intensity of $\eta_{0}=1.0$ starting
at $x_{0}=-\pi$.

Figure 6 illustrates the temporal evolution of the system governed by the
Advection-Diffusion equation. The initial condition is represented by the
dashed line profile. The advection term transports the system from left to
right, while the diffusion term smooths the initial distribution, causing it
to approach a Gaussian distribution more closely.

Figure 6: The time evolution of the diffusion-advection equation profile with
$D=0.01$ and $u=5$. The different colored solid lines represent
$t=0,0.1,0.2,\ldots,1.0$ with stepsize of $h=0.001$. The initial profile is a
top-hat function represented by the dashed line.

###  A.2 1D Burgers’ Equation

Burgers’ equation describes the dynamics of a fluid with nonlinear advection
and diffusion being commonly used as a model for a range of physical
phenomena, including turbulence, shock waves, and traffic flow. Here, we apply
our approach to solve the same problem as an example of the algorithm. The
Burgers’ equation is defined as

$\displaystyle\frac{\partial\eta}{\partial
t}=-\eta\frac{\partial\eta}{\partial x}+\nu\frac{\partial^{2}\eta}{\partial
x^{2}}$ (27)

with $\nu$ being the viscosity coefficient and $\eta$ being the velocity
field. The nonlinear term can be rewritten as
$\eta\partialderivative*{\eta}{x}=\tfrac{1}{2}\partialderivative*{\eta^{2}}{x}$
such that the pseudo-spectral formulation is defined by Eq. (11) with the
operators written as

$\displaystyle\widetilde{\mathcal{L}}_{{\bf\it k}}=-\nu
k^{2}\quad\text{and}\quad\widetilde{\mathcal{N}}_{{\bf\it
k}}(\eta)=-\frac{1}{2}ik\mathcal{F}\quantity{\eta^{2}}_{{\bf\it k}}.$ (28)

As an example, we set the viscosity coefficient as $\nu=0.001$. The size of
the system is $L=2\pi$ with the number of gridpoints $N=2^{12}=4096$. The
initial profile is a Gaussian given by $\eta_{0}(x)=\exp(-10(x/2)^{2})$.

Figure 7: The time evolution of the Burgers’ equation profile with
$\nu=0.001$. The different colored solid lines represent
$t=0,0.1,0.2,\ldots,1.0$ with stepsize of $h=0.001$. The initial profile is a
Gaussian function represented by the dashed line.

Figure 7 illustrates the time evolution of the $\eta$ field profile of the
Burgers’ equation. We can see the appearance of a shock front in the late
stages. As these shock waves propagate through the medium, they gradually
dissipate due to the viscosity effects.

##  Appendix B Python Codes

###  B.1 2D Cahn-Hilliard Model

⬇

1import numpy as np

2from scipy.fft import fft2, ifft2

3# Cahn-Hilliard model constants

4W = 1.0

5M = 1.0 # mobility

6kappa = 0.1 #gradient coeficient

7# Size of the system

8N = 2**8 # 2**8 = 256

9L = 16*np.pi

10x = np.linspace(0,L,N)

11dx = x[1]-x[0]

12# The time step definition

13h = 0.01

14T = 1500

15Nsteps = int(T/h)

16dframes = 1.0 # time step to output

17Nframes = int(T/dframes) #frames to the output

18nframes = Nsteps//Nframes

19# The array of outputs

20n = np.empty((Nframes,N,N), dtype=np.float32)

21# The Fourier variables

22n_k = np.empty((N,N), dtype=np.complex64)

23kx = np.fft.fftfreq(N, d=dx)*2*np.pi

24k = np.array(np.meshgrid(kx , kx ,indexing =’ij’), dtype=np.float32)

25k2 = np.sum(k*k,axis=0, dtype=np.float32)

26kmax_dealias = kx.max()*2.0/3.0 # The Nyquist mode

27# Dealising matrix

28dealias = np.array((np.abs(k[0]) < kmax_dealias )*(np.abs(k[1]) <
kmax_dealias ),dtype =bool)

29# The linear terms of PDE

30Loperator_k = -M*(kappa*k2**2+2*W*k2)

31# The non-linear terms of PDE

32def Noperator_func(n):

33 return -2*M*W*k2*fft2(-3*n**2+2*n**3)

34# Defining the time marching operators arrays

35# can be calculated once

36if method == ’IMEX’:

37 Tlinear_k = 1.0/(1.0-h*Loperator_k)

38 Tnon_k = dealias*h/(1.0-h*Loperator_k)

39elif method == ’IF’:

40 Tlinear_k = np.exp(h*Loperator_k)

41 Tnon_k = dealias*h*Tlinear_k

42elif method == ’ETD’:

43 Tlinear_k = np.exp(h*Loperator_k)

44 def myexp(x):

45 if x == 1: return 1.0

46 else: return (x-1.0)/np.log(x)

47 vmyexp = np.vectorize(myexp) # vectorize myexp (could be jitted)

48 Tnon_k = dealias*h*vmyexp(Tlinear_k)

49else: print(’ERROR: Undefined Integrator’)

50# Initial condition

51rng = np.random.default_rng(12345)

52noise = 0.02

53n0 = 0.5

54n[0] = n0 +noise*rng.standard_normal(n[0].shape)

55Noperator_k = n_k.copy() # auxiliary array

56nn = n[0].copy() # auxiliary array

57n_k[:] = fft2(n[0]) # FT initial condition

58# time evolution loop

59for i in range(1,Nsteps):

60 # calculate the nonlinear operator (with dealising)

61 Noperator_k[:] = Noperator_func(nn)

62 # updating in time

63 n_k[:] = n_k*Tlinear_k + Noperator_k*Tnon_k

64 # IFT to next step

65 nn[:] = ifft2(n_k).real

66 # test to output

67 if (i % nframes) == 0: n[i//nframes] = nn

###  B.2 2D Phase Field Crystal Model

⬇

1import numpy as np

2from scipy.fft import fft2, ifft2

3# PFC model constants

4r = -0.25

5M = 1.0 # mobility

6# Size of the system

7N = 2**8 # 2**8 = 256

8L = 16*np.pi

9x = np.linspace(0,L,N)

10dx = x[1]-x[0]

11# The time step definition

12h = 0.01

13T = 1500

14Nsteps = int(T/h)

15dframes = 1.0 # time step to output

16Nframes = int(T/dframes) #frames to the output

17nframes = Nsteps//Nframes

18# The array of outputs

19n = np.empty((Nframes,N,N), dtype=np.float32)

20# The Fourier variables

21n_k = np.empty((N,N), dtype=np.complex64)

22kx = np.fft.fftfreq(N, d=dx)*2*np.pi

23k = np.array(np.meshgrid(kx , kx ,indexing =’ij’), dtype=np.float32)

24k2 = np.sum(k*k,axis=0, dtype=np.float32)

25kmax_dealias = kx.max()*2.0/3.0 # The Nyquist mode

26# Dealising matrix

27dealias = np.array((np.abs(k[0]) < kmax_dealias )*(np.abs(k[1]) <
kmax_dealias ),dtype =bool)

28# The linear terms of PDE

29Loperator_k = -M*k2*(k2**2-2*k2+1+r)

30# The non-linear terms of PDE (with dealising)

31def Noperator_func(n):

32 return -(k2*M*fft2(n**3))

33# Defining the time marching operators arrays

34# can be calculated once

35if method == ’IMEX’:

36 Tlinear_k = 1.0/(1.0-h*Loperator_k)

37 Tnon_k = dealias*h/(1.0-h*Loperator_k)

38elif method == ’IF’:

39 Tlinear_k = np.exp(h*Loperator_k)

40 Tnon_k = dealias*h*Tlinear_k

41elif method == ’ETD’:

42 Tlinear_k = np.exp(h*Loperator_k)

43 def myexp(x):

44 if x == 1: return 1.0

45 else: return (x-1.0)/np.log(x)

46 vmyexp = np.vectorize(myexp) # vectorize myexp (could be jitted)

47 Tnon_k = dealias*h*vmyexp(Tlinear_k)

48else: print(’ERROR: Undefined Integrator’)

49# Initial condition

50rng = np.random.default_rng(12345)

51n0 = -0.085

52noise = 0.02*np.abs(n0)

53n[0] = n0 +noise*rng.standard_normal(n[0].shape)

54Noperator_k = n_k.copy() # auxiliary array

55nn = n[0].copy() # auxiliary array

56n_k[:] = fft2(n[0]) # FT initial condition

57# time evolution loop

58for i in range(1,Nsteps):

59 # calculate the nonlinear operator (with dealising)

60 Noperator_k[:] = Noperator_func(nn)

61 # updating in time

62 n_k[:] = n_k*Tlinear_k + Noperator_k*Tnon_k

63 # IFT to next step

64 nn[:] = ifft2(n_k).real

65 # test to output

66 if (i % nframes) == 0: n[i//nframes] = nn

## References

  * Chen [2002] L.-Q. Chen, Phase-Field Models for Microstructure Evolution, Annual Review of Materials Research 32, 113 (2002).
  * Gallego [2011] R. Gallego, Predictor-corrector pseudospectral methods for stochastic partial differential equations with additive white noise, Applied Mathematics and Computation 218, 3905 (2011).
  * Chen and Zhao [2022] L. Q. Chen and Y. Zhao, From classical thermodynamics to phase-field method, Progress in Materials Science 124, 100868 (2022).
  * Cahn and Hilliard [1958] J. W. Cahn and J. E. Hilliard, Free Energy of a Nonuniform System. I. Interfacial Free Energy, The Journal of Chemical Physics 28, 258 (1958).
  * Cahn [1959] J. W. Cahn, Free energy of a nonuniform system. II. Thermodynamic basis, The Journal of Chemical Physics 30, 1121 (1959).
  * Cahn and Hilliard [1959] J. W. Cahn and J. E. Hilliard, Free energy of a nonuniform system. III. Nucleation in a two-component incompressible fluid, The Journal of Chemical Physics 31, 688 (1959).
  * Elder and Grant [2004] K. R. Elder and M. Grant, Modeling elastic and plastic deformations in nonequilibrium processing using phase field crystals, Physical Review E 70, 051605 (2004), arXiv:0306681 [cond-mat] .
  * Elder  _et al._ [2007] K. R. Elder, N. Provatas, J. Berry, P. Stefanovic, and M. Grant, Phase-field crystal modeling and classical density functional theory of freezing, Physical Review B 75, 064107 (2007), arXiv:0607419 [cond-mat] .
  * Emmerich  _et al._ [2012] H. Emmerich, H. Löwen, R. Wittkowski, T. Gruhn, G. I. Tóth, G. Tegze, and L. Gránásy, Phase-field-crystal models for condensed matter dynamics on atomic length and diffusive time scales: an overview, Advances in Physics 61, 665 (2012).
  * Biner [2017] S. B. Biner,  _Programming phase-field modeling_ (Springer, 2017).
  * Boyd [2001] J. P. Boyd,  _Chebyshev and Fourier spectral methods_ (Courier Corporation, 2001).
  * Shizgal [2015] B. Shizgal, Spectral methods in chemistry and physics, Scientific Computation. Springer (2015).
  * Mortensen and Langtangen [2016] M. Mortensen and H. P. Langtangen, High performance Python for direct numerical simulations of turbulent flows, Computer Physics Communications 203, 53 (2016), arXiv:1602.03638 .
  * Giada  _et al._ [2002] L. Giada, A. Giacometti, and M. Rossi, Pseudospectral method for the Kardar-Parisi-Zhang equation, Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 65, 1 (2002), arXiv:0112333 [cond-mat] .
  * Ayub  _et al._ [2019] S. Ayub, H. Affan, and A. Shah, Comparison of operator splitting schemes for the numerical solution of the Allen-Cahn equation, AIP Advances 9, 125202 (2019).
  * Yoon  _et al._ [2020] S. Yoon, D. Jeong, C. Lee, H. Kim, S. Kim, H. G. Lee, and J. Kim, Fourier-Spectral Method for the Phase-Field Equations, Mathematics 8, 1385 (2020).
  * Chen and Shen [1998] L. Chen and J. Shen, Applications of semi-implicit Fourier-spectral method to phase field equations, Computer Physics Communications 108, 147 (1998).
  * Tegze  _et al._ [2009] G. Tegze, G. Bansel, G. I. Tóth, T. Pusztai, Z. Fan, and L. Gránásy, Advanced operator splitting-based semi-implicit spectral method to solve the binary phase-field crystal equations with variable coefficients, Journal of Computational Physics 228, 1612 (2009).
  * Martínez-Agustín  _et al._ [2022] F. Martínez-Agustín, S. Ruiz-Salgado, B. Zenteno-Mateo, E. Rubio, and M. Morales, 3D pattern formation from coupled Cahn-Hilliard and Swift-Hohenberg equations: Morphological phases transitions of polymers, bock and diblock copolymers, Computational Materials Science 210, 111431 (2022).
  * Milewski and Tabak [1999] P. A. Milewski and E. G. Tabak, Pseudospectral procedure for the solution of nonlinear wave equations with examples from free-surface flows, SIAM Journal of Scientific Computing 21, 1102 (1999).
  * Fornberg and Driscoll [1999] B. Fornberg and T. A. Driscoll, A Fast Spectral Algorithm for Nonlinear Wave Equations with Linear Dispersion, Journal of Computational Physics 155, 456 (1999).
  * Kassam and Trefethen [2005] A. K. Kassam and L. N. Trefethen, Fourth-order time-stepping for stiff PDEs, SIAM Journal of Scientific Computing 26, 1214 (2005).
  * Krogstad [2005] S. Krogstad, Generalized integrating factor methods for stiff PDEs, Journal of Computational Physics 203, 72 (2005).
  * Beylkin  _et al._ [1998] G. Beylkin, J. M. Keiser, and L. Vozovoi, A New Class of Time Discretization Schemes for the Solution of Nonlinear PDEs, Journal of Computational Physics 147, 362 (1998).
  * Cox and Matthews [2002] S. Cox and P. Matthews, Exponential Time Differencing for Stiff Systems, Journal of Computational Physics 176, 430 (2002).
  * Whalen  _et al._ [2015] P. Whalen, M. Brio, and J. V. Moloney, Exponential time-differencing with embedded Runge-Kutta adaptive step control, Journal of Computational Physics 280, 579 (2015).
  * Provatas and Elder [2011] N. Provatas and K. Elder,  _Phase-field methods in materials science and engineering_ (John Wiley & Sons, 2011).
  * Pope [1963] D. A. Pope, An exponential method of numerical integration of ordinary differential equations, Communications of the ACM 6, 491 (1963).
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Array programming with NumPy, Nature 585, 357 (2020).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods 17, 261 (2020).
  * Orszag [1971] S. A. Orszag, On the Elimination of Aliasing in Finite-Difference Schemes by Filtering High-Wavenumber Components, Journal of the Atmospheric Sciences 28, 1074 (1971).