# Accelerating solutions of one-dimensional unsteady PDEs with GPU-based swept
time-space decomposition

Daniel J. Magee  Kyle E. Niemeyer<EMAIL_ADDRESS>School of
Mechanical, Industrial, and Manufacturing Engineering, Oregon State
University, Corvallis, OR, USA

###### Abstract

The expedient design of precision components in aerospace and other high-tech
industries requires simulations of physical phenomena often described by
partial differential equations (PDEs) without exact solutions. Modern design
problems require simulations with a level of resolution difficult to achieve
in reasonable amounts of time—even in effectively parallelized solvers. Though
the scale of the problem relative to available computing power is the greatest
impediment to accelerating these applications, significant performance gains
can be achieved through careful attention to the details of memory
communication and access. The swept time-space decomposition rule reduces
communication between sub-domains by exhausting the domain of influence before
communicating boundary values. Here we present a GPU implementation of the
swept rule, which modifies the algorithm for improved performance on this
processing architecture by prioritizing use of private (shared) memory,
avoiding interblock communication, and overwriting unnecessary values. It
shows significant improvement in the execution time of finite-difference
solvers for one-dimensional unsteady PDEs, producing speedups of for a range
of problem sizes, respectively, compared with simple GPU versions and compared
with parallel CPU versions. However, for a more sophisticated one-dimensional
system of equations discretized with a second-order finite-volume scheme, the
swept rule performs worse than a standard implementation for all problem
sizes.

###### keywords:

GPU computing , partial differential equations , computational fluid dynamics
, high-performance computing , communication-avoiding algorithms , domain
decomposition

##  1 Introduction

High-fidelity computational fluid dynamics (CFD) simulations are essential for
developing aerospace technologies such as rocket launch vehicles and jet
engines. This project aims to accelerate such simulations to approach real-
time execution—simulation at the speed of nature—in accordance with the high-
performance computing development goals set out in the CFD Vision 2030 report
[1]. Classic approaches to domain decomposition for parallelized, explicit,
time-stepping partial differential equation (PDE) solutions incur substantial
computational performance costs from the communication between nodes required
every timestep. This communication cost consists of two parts: latency and
bandwidth, where latency is the fixed cost of each communication event and
bandwidth is the variable cost that depends on the amount of data transferred.
Latency in inter-node communication is a fundamental barrier to this goal, and
advancements to network latency have historically been slower than
improvements in other computing performance barriers such as bandwidth and
computational power [2]. Performance may be improved by avoiding external node
communication until exhausting the domain of dependence, allowing the
calculation to advance multiple timesteps while requiring a smaller number of
communication events. This idea is the basis of swept time-space decomposition
[3, 4].

Extreme-scale computing clusters have recently been used to solve the
compressible Navier–Stokes equations on over 1.97 million CPU cores [5]. The
monetary cost, power consumption, and size of such a cluster impedes the
realization of widespread peta- and exa-scale computing required for real-
time, high-fidelity, CFD simulations. While these are significant challenges,
they also provide an opportunity to develop new tools that increase the use of
the available hardware resources. As the authors of CFD Vision 2030 note,
“High Performance Computing (HPC) hardware is progressing rapidly and is on
the cusp of a paradigm shift in technology that may require a rethinking of
current CFD algorithms and software” [1]. Using graphics processing unit (GPU)
hardware as the primary computation device or accelerator in a heterogeneous
system is a viable, expedient option for high-performance cluster design that
helps mitigate these problems. For this reason, GPUs and other emerging
coprocessor architectures are increasingly used to accelerate CFD simulations
[6].

GPU technology has improved rapidly in recent years; in particular, NVIDIA
GPUs have progressed from Kepler to Pascal architecture in four years. This
development doubled and tripled peak single- and double-precision performance,
respectively [7]. In addition, the presence of GPUs in clusters, such as
ONRL’s Titan supercomputer, has become increasingly common in the last decade.
These advances have driven the development of software capable of efficiently
using and unifying the disparate architectures [8]. Although the ultimate
motivation of our work is accelerating the solution of PDEs—particularly
relevant to fluid flow—on distributed-memory systems, in this work we focused
on a single GPU-accelerated node/workstation in the design of the algorithms
and associated software. By investigating the effectiveness of the swept rule
on a workstation, we provide results that can be applied to simulations on a
single machine as well as an initial framework for understanding the
performance of the swept rule on heterogeneous computing systems.

The swept rule operates on a simple principle: do the most work possible on
the values closest to the processor before communicating. In practice, this
directive results in an algorithm that advances the solution in time at all
spatial points using locally accessible stencil values at the previous
timestep. Because the data closest to the processor is also the least widely
accessible, the strict application of this principle does not always provide
the best performance, but it is a useful heuristic for implementing the
procedure and analyzing its performance.

This study presents an investigation of the performance characteristics of
three swept rule implementations for a single GPU in a workstation. These
procedures are tested on three one-dimensional PDEs with numerical schemes of
varying complexity and compared with the performance of parallel CPU
algorithms and unsophisticated GPU versions. The (next) Section 2 describes
recent work on partitioning schemes for PDEs and communication-avoiding
algorithms, especially as applied to GPUs. Section 3 gives a brief overview of
the GPU architecture, particularly the thread and memory hierarchies. Section
4 discusses the swept rule, and our adjustments to the original algorithm in
response to the details of GPU architecture. Section 5 describes the swept
rule implementation in detail. In Section 6 we present the results of the
tests and, lastly, draw further conclusions in Section 7.

##  2 Related work

Alhubail and Wang introduced the swept rule for explicit, time-stepping,
numerical schemes applied to PDEs [9, 3, 4], and our work takes their results
and ideas as its starting point. The swept rule is closely related to cache
optimization techniques, in particular those that use geometry to organize
stencil update computation such as parallelograms [10] and diamonds [11]. The
diamond tiling method presented by Malas et al. [11] is similar to the swept
rule but uses the data dependency of the grid to improve cache usage rather
than avoid communication. Concepts such as stencil optimization using domain
decomposition on various architectures that are fundamental to this study are
explored by Datta et al. [12]. Their work explores comparisons between
parallel GPU and CPU architectures and tunes the stencil algorithm with nested
domain decomposition. The swept rule also has elements in common with
parallel-in-time and communication-avoiding algorithms.

Parallel-in-time methods [13], such as multigrid-reduction-in-time (MGRIT)
algorithms [14], accelerate PDE solutions with time integrators that overcome
the interdependence of solutions in the time domain, allowing parallelization
of the entire space-time grid. These methods calculate the solution over the
space-time domain using a coarse grid and iterate over successively finer
grids to achieve the desired accuracy. The use of coarse grids in parallel-in-
time methods reduces efficiency and accuracy when applied to nonlinear systems
[3]. This shortcoming is intuitive: since chaotic, nonlinear systems may
suddenly change in time, and coarse grids are prone to aliasing, the required
grid granularity diminishes gains in performance. The swept rule arises from
the same motivation, but does not seek to parallelize the computation in time
or vary dimensions during the process.

The swept rule does not alter the numerical scheme; it decomposes the domain
and organizes computation. That is, compared to a classic domain
decomposition, the swept rule performs the same operations in a different
order and location. In this way communication-avoiding algorithms share many
implementation details with the swept rule. Recent developments in
communication-avoiding algorithms for GPUs have generally focused on
applications involving matrices such as QR and LU factorization. The LU
factorization algorithm presented by Baboulin et al. [15] is motivated by the
increasing use of GPU accelerators in large-scale, heterogeneous clusters.
This method splits tasks between the GPU and CPU, minimizing communication
between devices. This allows the communication and computation performed on
each device to overlap, so all data transfer occurs asynchronously with
computation. We explore this approach—overlapping data transfer with hybrid
computation—in this article. The motivation and structure of our study is
comparable to the work of Anderson et al [16]. They developed methods for
arranging and tuning computation for single general-purpose GPU (GPGPU) in a
desktop workstation, without altering the basic QR factorization algorithm.
Similarly, this study focuses on adapting the swept rule to a single GPU. The
swept rule is a strategy for arranging the computational path of explicit
numerical methods, and this work seeks to design the data structures and
operations used in that path to achieve the best performance on GPUs.

##  3 GPU architecture and memory

The impressive parallel processing capabilities of modern GPUs resulted from
architecture originally designed to improve visual output from a computer. GPU
manufacturers and third-party authors [17, 18, 19, 20] have described the
features of this architecture in great detail, while others discussed general
best practices for developing efficient GPU-based algorithms [21, 6].
Particular aspects of this architecture, such as the unique and accessible
memory hierarchy, are at the core of this work, so some explanation of its
relevant elements is necessary before describing the details of the
implementation.

Programs that run on the GPU can be implemented using several software
packages, the most common of which are the OpenCL and OpenACC frameworks, and
the CUDA parallel computing platform. These packages use different
nomenclatures and are compatible with different hardware types. In this
project all programs use CUDA, which is exclusively compatible with NVIDIA
GPUs; therefore, all descriptions of GPU hardware presented here use the CUDA
nomenclature. CUDA programs consist of functions, referred to as kernels,
launched from a C/C++ host program. The CPU executes the host code and
specifies the size and number of blocks when calling a kernel, the stream
(queue) in which the kernel will be launched, and the amount of shared memory
to be allocated per block at runtime. All threads in a warp—a group of $32$
threads that execute as a single-instruction multiple thread (SIMT) unit—must
be in the same block, so for good practice blocks should launch with some
multiple of 32 threads.

The information presented here is valid for all NVIDIA GPUs with compute
capability $3.0$ or higher (i.e., Kepler architecture or later). The device
used in this study is a Tesla K40c GPGPU, compute capability 3.5. This device
contains 15 streaming multiprocessors, each capable of processing 64 warps of
32 threads, or 2048 total threads, at once [7]. A maximum of $16$ blocks may
concurrently reside on a streaming multiprocessor; blocks may not be split
between streaming multiprocessors. While each streaming multiprocessor can
support 2048 resident threads, in practice their capacity is often lower
because each thread or block makes demands on limited memory resources—most
notably the shared memory and registers. Each streaming multiprocessor on the
Tesla K40c has $48\text{\,}\mathrm{kB}\text{/}$ of shared memory and 65536
registers available. Registers offer the fastest access, but are the most
limited memory type and are private to each thread, but can be accessed by
other threads in the same warp using shuffle operations available on devices
with compute capability 3.5 or higher. Shared memory is a controllable portion
of the L1 cache accessible only to threads within a block. As a result, for a
thread to read a value stored in shared memory in a different block, a thread
with access to that value must write the value to global memory where the
reader thread has access. Global memory is the slowest and most plentiful
memory type, and where data copied from the host program resides. Global
memory stores all variables passed to a kernel and large arrays declared
therein.

Other memory types in the CUDA memory hierarchy include constant, texture, and
surface; of these, the work presented here only uses constant memory. Constant
memory is read-only, available to all kernels for the lifetime of an
application, and quick to access when all threads access the same location.
This makes it a convenient and performance conscious choice for storing
constant values of the governing equations calculated at runtime [18].

##  4 Methodology

###  4.1 Experimental method

The primary goal of this study is to compare the performance of the swept rule
to a simple domain decomposition scheme, referred to as Classic, on a GPU. A
domain decomposition scheme is a way of splitting up a large problem so tasks
can be shared by discrete workers working on distinct parts of the problem. In
this case the decomposition scheme divides the space-time domain of the PDE
solution. We will compare the performance of these methods by the time cost of
each timestep, including the time required to transfer data to/from the GPU.
While encoding the Classic is relatively straightforward, finding the best
approach for the swept rule on the GPU presents a more subtle problem.

In the original swept rule approach [3], the spatial domain is partitioned
into independent pieces called “nodes” that correspond to compute nodes on a
distributed system with private memory spaces. A major concern in adapting
this nodal analogy to a single GPU is the type of memory allocated for the
working array, the container for the values that are the solution to and the
basis for each timestep. Several available approaches exist to map the
original analogy for a node to a single GPU; here, we will explore three of
them: Shared, Hybrid, and Register (which we will describe in detail in
Section 5).

In all approaches we map one thread to one spatial point. Handling more than
one spatial point per thread would allow for larger nodes but would require
more resources and complicate the procedure without reducing thread idleness.
Additional swept rule properties could be adjusted for potential
implementation variants such as the data structure to hold the working values,
and the method to globally synchronize threads. However, for the purposes of
this study, we did not vary these attributes. In all cases in this study the
working array is a standard, one-dimensional C array with two flattened rows;
and kernel calls are implicitly synchronized by returning control to the queue
in the host program.

###  4.2 First-order domain of dependence

(a) The first step of the swept rule. Values at $t=0$ are split between nodes
$0$ and $1$, which compute solutions in their domain of dependence, a triangle
in the space-time plane. The edge values are collected in global arrays
$L_{0}$/$R_{0}$ and $L_{1}$/$R_{1}$.

(b) The second step in the swept rule. The nodes pass their right edge to the
neighboring node. The passed values become the initial left edge, and the left
edge from the previous stage becomes the right edge. Each node advances
through their domain of dependence, a diamond in space-time.

Figure 1: The first two steps of the swept rule for a numerical scheme with a
first-order domain of dependence. $L_{i}/R_{i}$ refer to left/right arrays of
node $i$, which collect edge values shown with thick bordered dots [22].

A domain of dependence is a region on the space-time grid that can be
completed by a numerical scheme with some set of initial values. For the
purposes of this project, numerical schemes consist of stencil operations,
where a value at a grid point is updated based on the weighted contribution of
values at grid points in the vicinity. A three-point stencil uses values at
neighboring grid points only, and the timestep of any numerical scheme can be
decomposed into a series of sub-timesteps that only require a three-point
stencil [23]. The numerical scheme defines the order of the domain of
dependence: it increases by one for every two sub-timesteps required per
timestep. Therefore, a domain of dependence is _first-order_ if all sub-
timesteps in the numerical scheme use a three-point stencil, and intermediate
values are required no more than two steps after they are calculated. The
initial incarnation of the swept rule presented by Alhubail and Wang [3, 9]
decomposes multi-step timesteps and large stencils into sub-timesteps with a
three-point stencil. This regularizes the procedure and ensures that all
equations and schemes can be evaluated using a swept decomposition with a
first-order domain of dependence, but requires more memory to store the
intermediate values that result from each sub-timestep. In order to conserve
limited private memory, for the GPU-based swept rule a first-order domain of
dependence is applicable to schemes that require two or fewer three-point
stencil sub-timesteps per timestep. Figure 1 shows the first two stages of the
swept rule using $k=2$ nodes with $n=16$ spatial points.

At the start of the first step, shown in Figure 1(a), the initial conditions
are passed to the kernel and each node evaluates the solution at as many
points as possible in the space-time grid. The initial domain of dependence
forms a triangle. The working array stores the solutions as each node steps
through time and is maintained in a fast memory space private to node member
threads. When each node cannot advance any further it’s necessary for each to
pass one edge to a neighboring node and retain the values of the other edge.
Figure 1(b) shows how the second step proceeds from the first using the edge
values.

(a) At the end of a kernel (completion of an inverted triangle), the working
array is stored and passed from private memory to global memory. The location
shows the index of the working array in (private) shared memory. The number in
the shape refers to the offset global memory index where the working value is
passed.

(b) Edges are reinserted to seed the inverted triangle of a new swept cycle.
The location shows the index of the working array in (private) shared memory
receiving the global memory value denoted by the number.

Figure 2: The procedure for edge passing shown in Figure 1(b). The global
arrays $L_{i}$ and $R_{i}$ are represented by circles and squares,
respectively, and the numbers in those shapes represent the position of the
value in the global array. The axes describe the location of the value in the
working array [22].

For the nodes to communicate as Figure 1 illustrates, the values on the edges
must be passed to global arrays available to all nodes since any sufficiently
fast memory location is private to a subset of threads terminated on kernel
exit. Consequently, any data necessary to continue the computation must be
shared between nodes using an intermediary container. The information needed
to begin a new nodal cycle, including the values retained by each node, must
be read in from, and out to, global memory at the beginning and end of each
kernel, respectively. Figure 2(a) shows how the working array values are
stored and passed to the global arrays.

Figure 1(b) illustrates the reason the two edges are stored individually as
$L_{i}$ and $R_{i}$: only one of the edges is passed between nodes. After the
first step, the right edge is passed to its right neighbor node; the left edge
is stationary. At the beginning of the subsequent kernel, $L_{i}$ and $R_{i}$
are swapped and reinserted into the working array to seed the next progression
with one data transfer event from global to private memory as shown in Figure
2(b). When the right edge is passed between nodes, node $0$ is split across
the spatial boundary and must apply the boundary conditions at its center.

The swept rule allows the computation to advance by $n$ timesteps with two,
rather than $n$, global memory accesses, where $n$ is the number of threads
per block (or the number of spatial points per node). The procedure advances
by passing values in the alternating directions and zig-zagging the location
of the nodes in this fashion until the simulation is complete. Since the
diamonds shown in Figure 1(b) do not store all the values at a single
timestep, the simulation can only output values when a complementary triangle
is computed and the final $n$-length local tier is returned. This kernel can
only be called after the values are passed to the left, so the results can
only be read out every $n$th timestep.

Although we already described the working array and showed in Figure 2 how the
relevant values are communicated between the nodes, it is instructive to
outline the performance concerns that motivate this arrangement. As Section 3
describes, the number of resident threads on a streaming multiprocessor
depends on the GPU architecture and resources requested at kernel launch. For
instance, storing every double-precision value in the triangle shown in Figure
1(a) in shared memory in a kernel with $512$ threads per block would require
$8\times 65792=$514\text{\,}\mathrm{kB}\text{/}$$ of shared memory—this is
over 100 times greater than the $48\text{\,}\mathrm{kB}\text{/}$ limit for
NVIDIA GPUs with Kepler architecture. The maximum number of threads per
streaming multiprocessor would be limited to 128, which would negatively
impact program performance.

Figure 1(a) shows that the interior of the triangle is only needed to progress
to the next timestep, and that edges on even and odd tiers do not overlap in
the spatial domain. Thus, the triangle may be stored as a matrix with two
rows, where the first and second rows contain the even and odd sub-timesteps
results, respectively. The interior values are overwritten once they are used,
and only the edge values remain. Figure 2(a) shows the result of a local
computation with this method. The last two values, the tips of each triangle,
are copied into both arrays.

At the start of the next kernel, the left and right arrays are inserted into
the working array to seed successive calculations as Figure 2(b) shows. The
edges of the previous cycle’s first row are moved to the center, and the
center of the current top row is now left open for the first tier of the
inverted triangle. Each row requires space for $n+2$ values because two edge
values are required on either side to complete the stencil for the longest row
where $n$ values are computed. The computation proceeds by filling the empty
two indices on the top row, overwriting the bottom row’s middle four indices,
and so on. In contrast to the memory demands of storing the entire nodal
computation, this method uses only $2\times(n+2)$ values. For a block with n =
512 threads, this requires $8\times
1028\approx$8\text{\,}\mathrm{kB}\text{/}$$. In general, a stencil with width
$k$ would require storing $2\times(n+k-1)$ values in the working array. By
reducing the amount of shared memory to only $8\text{\,}\mathrm{kB}\text{/}$,
the kernel is less limited and achieves higher occupancy, the number of
threads each streaming multiprocessor is capable of handling simultaneously.

###  4.3 Higher-order domain of dependence

The first-order domain of dependence suffices for relatively simple problems.
However, more complicated problems with elements such as nonlinear equations,
discontinuities, or higher-order derivatives require more sophisticated
procedures. The original swept rule program calculates and stores intermediate
values at sub-timesteps for higher-order schemes to avoid using larger
stencils [3]. Breaking a timestep into a series of sub-timesteps allows any
numerical scheme for any equation of the same dimension to be decomposed in
the same way since all stages in the computation depend on the minimum
stencil, that is, only the neighbors of the current spatial point [23]. For
example, a second-order in time midpoint method applied to a fourth-order
differential equation would require four sub-timesteps per timestep. The first
sub-timestep would find the second derivative of the dependent variable so
that the second sub-timestep, the midpoint solution, would only require a
three point stencil: the second derivative and initial values at the
neighboring spatial points. The third sub-timestep would find the second
derivative using the midpoint solution, and the final sub-timestep would
complete the timestep using the results of the second and third sub-timesteps
on the three-point stencil and the previous timestep solution at the current
spatial point. This approach presents a storage and data transfer problem on
the GPU because values in the interior of the working array are overwritten
two sub-timesteps after they are calculated. These forgotten but required
values are marked with an “$\times$” in Figure 3(b).

(a) First step of swept rule for discretization with four, three-point stencil
sub-timesteps per timestep.

(b) Using four sub-timesteps per timestep will overwrite values marked with an
“$\times$” before they are needed.

Figure 3: Conflicts in the domain of dependence for discretizations requiring
more than two sub-timesteps per timestep [22].

Saving four values per tier would fix the problem, but requires a larger
matrix in shared memory, which would diminish occupancy and require more
unnecessary values to be passed between nodes. Figure 4 shows our solution to
this problem: a five-point stencil that requires two sub-timesteps per
timestep—a predictor and a final step. This flattens the triangle or diamond
in the time domain and requires four values per sub-timestep, but the two-row
matrix may be used as described in Section 4.2 and illustrated by Figure 2
with minor adjustments. The same number of values are transferred between
nodes in each communication, but inevitably more communication events are
required to advance the solution. Conveniently, the predictor-corrector method
ensures that all odd tiers, the second matrix row, will contain predictor
values, and the bottom row will hold final values.

(a) First step in swept rule for second-order domain of dependence. The
working array is able to be folded and passed as shown in Figure 2

(b) Second step in the swept rule second-order domain of dependence.

Figure 4: The first two steps of the swept rule for a numerical scheme with a
second-order domain of dependence [22].

The problem that motivates our adjustment to the swept rule is the result of
both the midpoint method and the five-point stencil discretization. Either of
these circumstances alone would accommodate the method described in the
previous section. Generally, the decomposition must be adjusted in this
fashion when there are more than two sub-timesteps per timestep. It would need
to be adjusted further for problems with more than four sub-timesteps per
timestep.

##  5 Implementation

###  5.1 Swept rule variants

While the order of the swept domain of dependence is a natural consequence of
the equations and numerical scheme, the computational hierarchy available on
GPUs allows this structure to be implemented in different ways. To more
thoroughly investigate the effects of the swept rule, we consider three
versions: Shared, Hybrid, and Register, and compare them with a basic
decomposition algorithm, Classic. In all of these programs, one GPU thread is
assigned to one spatial point.

⬇

__global__ void

upTriangle(const REAL *IC, REAL *right, REAL *left) {

//tid = threadIdx.x (bottom row idx),

//tidT = tid+blockDim.x (top row idx);

//tids and tidTs: stencil for each tidT and tid respectively.

//shareT = working array in shared memory

//Read initial data in from IC to shareT.

if (tid > 1 && tid <(blockDim.x-2)) {

shareT[tidT]=predictorStep(shareT, tids);

}

__syncthreads();

//4 is stencil length - 1

for (int k = 4; k<(blockDim.x/2); k+=4)

{

if (tid < (blockDim.x-k) && tid >= k) {

shareT[tid]+=finalStep(shareT, tidTs);

}

k += 2;

__syncthreads();

if(tid < (blockDim.x-k) && tid >= k) {

shareT[tidT]=predictorStep(shareT,tids);

}

__syncthreads();

}

//Read out the edges to left and right arrays.

}

Figure 5: Main loop of the starting kernel for the swept rule as illustrated
by Figure 4(a).

⬇

__global__ void

classicKS(const REAL *ks_in, REAL *ks_out, bool final)

{

//Global Thread ID

int gid = blockDim.x * blockIdx.x + threadIdx.x;

//number of spatial points - 1

int lastidx = ((blockDim.x*gridDim.x)-1);

//Stencil indices.

int gids[5];

//True for all spatial points from periodic BCs.

#pragma unroll

for (int k = -2; k<3; k++) {

gids[k+2] = (gid + k) & lastidx;

}

//Final is false for predictor step, true otherwise.

if (final) {

ks_out[gid] += finalStep(ks_in, gids);

}

else {

ks_out[gid] = predictorStep(ks_in, gids);

}

}

Figure 6: Classic kernel for Kuramoto-Sivashinsky solver. Final and predictor
step functions can be written in C with __device__ keyword.

The Classic algorithm is a naive GPU implementation of the numerical solution
and the baseline against which the efficacy of the swept rule is measured. It
advances one sub-timestep per kernel call and uses global memory to store the
working array. Figure 6 shows the structure of the procedure, and it is
similar for all problems and discretizations.

We consider the Shared strategy for implementing the swept rule the most
natural way to map the analogy of CPU nodes to GPU architecture. It is applied
to every test case and is considered the “default” swept rule GPU version for
comparing the performance of the GPU programs with their MPI-based CPU
counterparts [9]. The Shared version treats a block as a node and uses shared
memory for the working array. Each block has exclusive access to a shared
memory space and may contain up to $1024$ threads, so it has access to fast
memory and the capacity for various node sizes. There are some drawbacks to
this version: a high number of idle threads for sub-timesteps where the domain
of dependence contains relatively few spatial points, poor utilization of CPU
resources, and the fact that shared memory is not the fastest memory type
[24].

The Hybrid strategy uses the same GPU procedure as Shared, uses the CPU to
compute the node that is split across the boundary, as seen in Figure 1(b),
Transfers between the host and device are costly operations, but the devices
can execute instructions concurrently—so if the CPU can complete the boundary
node before the GPU finishes the other nodes, no penalty arises. In this
study, we apply this strategy to problems with non-periodic boundary
conditions as a way to mitigate the underutilization of the CPU and the thread
divergence that results from applying boundary conditions in a GPU kernel.

The Register approach is applied to problems with periodic boundary conditions
because of the difficulty involved in applying boundary conditions using warp
shuffle functions. This implementation limits the number of points in a node
to the size of a warp, $32$ threads (which has been constant over several
iterations of NVIDIA GPUs). In this version the values are initially read from
global memory to shared memory to the registers. Passing the values to the
intermediate shared memory is necessary because the shuffle operations that
trade registers between threads only operate on active threads in the warp
being called; if some threads are masked, they will be unable to supply the
necessary stencil values. Thus, data must be moved between memory levels at
each tier rather than once at the start and end of the kernel. This seriously
limits the Register approach, but it still warrants exploration since
registers are the fastest memory type.

###  5.2 Test cases

We present three test cases to demonstrate the performance and functionality
of the GPU-based swept rule in one spatial dimension: the heat equation,
Kuramoto–Sivashinsky (KS) equation, and Euler equations for compressible flow.
Appendices Appendix B–Appendix D contain the full derivations of the
procedures for the numerical solutions. First, we chose the heat equation for
its simplicity and familiarity. Here it is discretized with a first-order
scheme using forward differencing in time and central in space. Next, we
selected the KS equation to demonstrate the swept rule for higher-order,
nonlinear PDEs. We discretized the KS equation with second-order, central
differencing in space, which requires a five-point stencil, and a second-order
Runge–Kutta method in time. Lastly we chose to solve the Euler equations, a
system of quasilinear, hyperbolic equations for describing compressible,
inviscid flow. The conservative form of these equations are applied to the Sod
shock tube problem to demonstrate the application of the swept rule to a
canonical CFD problem involving discontinuities and several dependent
variables. These equations are discretized with a second-order, finite-volume
scheme in space and a second-order method in time. The heat equation requires
a first-order domain of dependence, while the KS and Euler equations require
second-order domains of dependence. These problems also provide examples of
various types of boundary conditions. The heat and Euler equations are solved
with reflective and Dirichlet boundary conditions, respectively; therefore,
the boundary conditions must be imposed with control flow. The KS equation
uses periodic boundary conditions, which is a convenient formulation for the
swept rule that splits a node across the boundary.

##  6 Results and discussion

All tests presented here were performed on a single workstation with a Tesla
K40c GPU and an Intel Xeon 2630-E5 CPU with eight cores and 16 potential
threads. The GPU-based swept rule algorithms and test cases were implemented
1DSweptCUDA v2 [25]. For the results we present here, each program was
executed in double precision with $\\{2^{x}\mid x\in\mathbb{N}\mid 4<x<11\\}$
threads per block and $\\{2^{x}\mid x\in\mathbb{N}\mid 10<x<21\\}$ spatial
points. Each run advanced $50,000$ timesteps and recorded the average time per
timestep; the initial GPU memory allocations and data transfer between the
host and device are included in the overall time measurement. Then, we
collected the best time per timestep for each number of spatial points. We
repeated this procedure five times and took the average to obtain the results
presented here. There were no significant differences in the results between
the tests for the same configuration.

Figures 7(a) and 8(a) show the execution time per timestep for all algorithms
applied to the heat and KS equations. When applied to these problems, the
swept rule algorithm outperforms the Classic procedure, and the GPU-only
shared memory version, Shared, is faster than the alternate swept rule
versions. Figures 7(b) and 8(b) show the speedup of the swept rule programs,
or the ratio of time costs with the Classic version. Both cases exhibit
similar performance patterns: Shared generally provides a larger speedup for
small spatial domains ($<${10}^{4}$$ spatial points), but only a $2$$\times$
speedup for large ones ($>10^{5}$ spatial points). Figures 7(a) and 8(a) show
the performance trends of the algorithms with respect to the spatial domain
size. Their time costs are insensitive to the spatial domain at smaller domain
sizes and grow linearly with increasing domain size after about
$3\text{\times}{10}^{5}$ spatial points. The similarity in the performance
trends of the heat and KS programs is intuitive; although the KS equation is
nonlinear, fourth-order, and discretized with a higher-order scheme, both are
continuous, scalar equations and use finite difference schemes.

(a) Time cost of GPU classic and swept domain decomposition algorithms.

(b) Speedup of swept rule programs with respect to the Classic version.

Figure 7: Performance comparison of the GPU heat equation programs [22].

(a) Time cost of GPU classic and swept domain decomposition algorithms.

(b) Speedup of swept rule programs with respect to the Classic version.

Figure 8: Performance comparison of the GPU Kuramoto–Sivashinsky equation
programs [22].

The speedup of the swept rule declines as the number of spatial points in the
domain increases, caused primarily by occupancy and the fixed cost of kernel
launch. The occupancy of a kernel is the number of threads that can reside
concurrently on a streaming multiprocessor at launch. This quantity is
determined by the block size, since blocks of threads are indivisible, and the
resources requested by the kernel. The swept rule requests more resources,
specifically shared memory and registers, than Classic to store the working
array and carry out its procedure, which involves more steps. If more
resources are requested than are available to the streaming multiprocessors on
the device, the GPU launches waves of blocks. If the occupancy is limited by
resource allocation, the kernel must begin launching waves of blocks on
domains with fewer spatial points. Each wave must wait until the previous one
completes before beginning; theoretically this implies two waves will cost
about twice as much as one. As the number of spatial points in the domain
grows, more blocks are launched, and the difference in the number of waves per
kernel call increases more quickly for the swept rule program. This conclusion
arises from the observation that the time cost of the swept rule versions
begins growing linearly at a smaller spatial domain size, about
$2\text{\times}{10}^{4}$ points, than Classic, about $6\text{\times}{10}^{4}$
points, as shown in Figures 7(a) and 8(a).

By timing the launch of many empty kernels and taking their average, we
measured the fixed cost of kernel launch on the testing workstation to be
about $4\text{\,}\mathrm{\SIUnitSymbolMicro s}\text{/}$. We conclude that this
cost dominates the performance of Classic at small problem sizes because it is
quite close to the cost of each timestep for spatial domains with less than
$3\text{\times}{10}^{4}$ spatial points for both the heat and KS equations.
This fixed cost accounts for less time cost at larger spatial domain sizes
where all kernels must launch several waves of blocks, so the portion of the
swept rule speedup from avoiding kernel launches becomes negligible, and the
overall speedup of the swept rule is reduced.

Each of the swept rule variants experiences these trends, but in all cases the
Shared version performs better. Figure 7 shows the heat equation test, where
the Hybrid version performs as well as Shared for large spatial domain sizes
and 2–3 times worse for small ones. The Hybrid version uses the same GPU
kernels as Shared, but uses the CPU to calculate the first node when it is
split across the boundary. At smaller domain sizes, the data transfer between
the host and device dominates the cost of the Shared scheme. At larger domain
sizes, this cost drops in comparison with the overall costs, but at the same
time the thread divergence caused by handling boundary conditions also drops
in importance. Thus, the Shared and Hybrid methods perform similarly at larger
domain sizes (i.e., above $6\text{\times}{10}^{4}$ spatial points).

Figure 8 shows the KS equation test, where, similar to the heat equation, the
Shared version performs better in all cases. While registers are the fastest
memory type, because of memory access rules the Register version must use a
warp as a node, which limits the domain of dependence to 32 spatial points. A
Register program with $n=32$ must launch four times as many kernels to advance
the same number of timesteps as Shared with $n=128$ (most often the best block
size). Despite these differences, the Register version exhibits a similar
performance trend to Shared, which shows about $1.5$$\times$ the speedup of
Register at all spatial domain sizes. The limit on node size results in a
constant four timesteps per cycle for the KS equation, which only reduces the
number of communication events by $1/8$ compared with Classic.

(a) Time cost of GPU classic and swept domain decomposition algorithms. Hybrid
and Shared lines overlap.

(b) Speedup of swept rule programs with respect to the Classic version.

Figure 9: Performance comparison of the GPU Euler equation programs [22].

In contrast to the previous test cases, Figure 9 shows that the classic
decomposition outperforms the swept rule at all problem sizes when applied to
the Euler equations. This result differs from the findings of Alhubail and
Wang [3] for the swept rule implemented only with MPI on CPUs, which produces
roughly equivalent speedups for both the KS and Euler equations. The GPU
implementation of the swept rule for the Euler problem involves much greater
arithmetic intensity than the other problems, causing greater low-level memory
usage. This limits the performance of the swept scheme on a GPU more than a
CPU, reducing the benefits of the scheme over the classic approach within a
single GPU compared with the other problems. Those cases involve higher ratios
of floating-point operations to memory accesses (both read and write). This
will particularly degrade performance as the intra-domain timesteps proceed
and nodes become inactive.

With regard to the potential performance of the swept rule, this result is not
encouraging for more complex problems in more dimensions. However, on a larger
scale in a cluster, where global memory is used for the working array and the
communication to be avoided is through a physical network between processors,
we expect that the swept rule will provide a benefit for these types of
problems.

(a) Time cost of CPU and GPU programs for KS equation.

(b) Performance improvement of GPU program compared to the same algorithm on
the CPU.

Figure 10: Performance comparison of CPU (MPI) and GPU (CUDA) programs for the
KS equation [22].

Figure 10 compares the CPU-based KS equation programs with the GPU algorithms,
Classic and Shared. This CPU version is parallelized with MPI, similar to that
originally presented by Alhubail and Wang [9]. Figure 10 refers to the Shared
program as SweptGPU because it serves as the “default” swept rule version on
the GPU for comparison with the CPU-based swept rule. The original CPU program
was designed for a distributed-memory cluster, but the Alhubail and Wang’s
original performance study used only two processes on two CPUs [3]. The
parallel CPU program was tested similarly to the GPU programs: it was run on
the same workstation with the same spatial domain sizes and $216$ threads.
Here we compare the best results for each number of spatial points, which
usually occurred with 16 threads. This did not degrade the performance
compared with the original study; in fact, both CPU versions of the program
performed significantly better, and for most spatial domain sizes the CPU-
based swept rule improved by about . Since the test used up to eight times the
number of threads, this result supports the validity of repurposing this code
and comparing the result with the GPU versions.

Figure 10 illustrates the performance benefits of the GPU architecture more
than the swept rule itself. On small spatial domains (e.g., less than
$1\text{\times}{10}^{4}$ points) the GPU can assign one thread to each spatial
point and process all in a single wave. There it improves performance over the
CPU version by less than , because the work does not fully utilize the GPU. On
larger spatial domains (e.g., more than $1\text{\times}{10}^{4}$ points),
where the work completely utilizes the resources of the GPU, the performance
increases continue growing with problem size. The increase of the GPU
programs’ improvement with respect to the number of points in the spatial
domain inverts the trend of the swept rule’s improvement with respect to the
Classic kernel. Both GPU algorithm types show a speedup of about for the
smallest spatial domain. At the other end, the speedup grows to about for the
swept and for the classic domain decomposition.

##  7 Conclusions

In this study we compared the time per timestep of three swept rule time-space
decomposition implementations to a simple domain decomposition scheme,
Classic, for GPU accelerators. These programs were evaluated on a range of
spatial domain sizes. The Classic and Shared programs were also compared with
their CPU-based counterparts parallelized with MPI.

Generally, the swept rule performs better relative to Classic when the kernel
launch cost is a substantial element of the program cost and the spatial grid
has fewer points than concurrently launchable threads on the GPU. Once this is
exceeded, the launch cost penalty becomes negligible and greater resource
usage penalizes the swept rule in the form of reduced occupancy and reduced
availability of the L1 cache, which is reserved for shared memory. But, like
the initial performance penalty for the Classic program, the resource usage
penalty does not scale and all programs see their time cost rise nearly
linearly with respect to spatial domain size after about ${10}^{5}$ spatial
points.

Both alternate swept-rule procedures, Hybrid and Register, performed slightly
worse than the Shared version. We attribute the failure of the Hybrid routine
to improve performance to improvements made in the handling of the boundary
conditions in Shared rather than the cost of host-to-device communication,
which is handled with asynchronous streams. Though the Hybrid version does not
improve performance, it does not substantially degrade it either, especially
with large spatial domain sizes. Ultimately, Hybrid computation seeks to solve
a problem more efficient to solve within the GPU paradigm. The Register
computation shows more promise but falls short because of the limited node
size.

This study also shows that implementing a Classic decomposition on the GPU for
an explicit numerical scheme is simple and may result in noticeable
performance improvements. This may be enough for many applications, but in
performance-critical cases the swept rule may further reduce execution time.
The performance of the Classic decomposition programs may be improved by
altering the synchronization method. Implicit synchronization incurs the cost
of kernel launch at each sub-timestep, and this kernel executes so quickly at
small problem sizes that synchronization cost dominates the performance of the
program. Using off-the-shelf, GPU-based synchronization [26] could also
provide performance benefits for the swept rule, and in particular Register,
which is also limited by kernel launch cost for small spatial domains.

The specific results presented here depend on the hardware used. A commercial
GeForce GPU with Kepler architecture would perform worse than the Tesla K40c,
which is designed for computation as opposed to actual graphics. Despite the
significant performance increase shown here, the Tesla K40c is three-year-old
technology. The current state-of-the-art GPGPU with Pascal architecture, the
Tesla P100, offers twice the K40c base clock speed, nearly four times the
number of streaming multiprocessors, and an extra
$16\text{\,}\mathrm{kB}\text{/}$ of shared memory independent of the L1 cache.
Additional dedicated shared memory could dramatically impact the swept rule’s
performance for the Euler equations. We predict that this device could at
least halve the execution times shown here and maintain insensitivity to
problem size up to over ${10}^{5}$ spatial points, potentially resulting in
speedups on the order of ${10}^{3}$ over CPU parallel versions. Future GPU
accelerators could further improve performance, as well as devices with
similar SIMT architectures like Intel MIC/Xeon Phi.

Our future work will focus on further developing the swept rule for use on
distributed memory systems with heterogeneous node architectures and
implementing the swept rule in higher dimensions on the GPU. We expect the
performance benefits to increase from the swept rule when applied to a
distributed memory system comprised of nodes with CPUs and GPUs, where network
communications incur more latency than intra-GPU memory access. However, we
anticipate the performance benefits of the swept rule will diminish for higher
dimensions; in two dimensions the swept rule requires three stages of inter-
node communication to advance one cycle, analogous to a single diamond in one
dimension. This, along with the associated increase in required memory
allocation, arithmetic intensity, and kernel launch events could limit the
performance of the swept rule as observed by Alhubail et al. [4]. In spite of
these challenges, recent developments such as grid-wide synchronization and
increased shared memory capacity could aid the performance of the swept rule
and permit the exploration of additional design options.

##  Appendix Appendix A Availability of materials

The results for this article were obtained using 1DSweptCUDA v2 [25]. All
figures, and the data and plotting scripts necessary to reproduce them, are
available openly under the CC-BY license [22].

##  Appendix Appendix B Heat equation

The unsteady heat conduction equation without volumetric heat flux is

$\frac{\partial T}{\partial t}=\alpha\nabla^{2}T\;.$ (1)

where $T$ is temperature, $t$ is time, and $\alpha$ is thermal diffusivity. In
one dimension, this is reduced to

$\frac{\partial T}{\partial t}=\alpha\frac{\partial^{2}T}{\partial x^{2}}\;.$
(2)

Discretizing Eq. (2) with forward differencing in time and central in space
yields

$\frac{T_{i}^{m+1}-T_{i}^{m}}{\Delta
t}=\alpha\frac{T_{i+1}^{m}+T_{i-1}^{m}+2T_{i}^{m}}{\Delta x^{2}}\;,$ (3)

where $i$ is the spatial node index and $m$ is the time index corresponding to
time $t^{m}$. This is a first-order, explicit, finite-difference
approximation. To step forward in time, define the Fourier number,
$\text{Fo}=\frac{\alpha\Delta t}{\Delta x^{2}}$, where $\Delta t$ is the
timestep size and $\Delta x$ is the spatial grid size, and solve for
temperature at the next timestep:

$T_{i}^{m+1}=\text{Fo}(T_{i+1}^{m}+T_{i-1}^{m})+(1-2\text{Fo})T_{i}^{m}\;.$
(4)

In this study the approximation is evaluated with insulated boundary
conditions at both ends and $n$ spatial points:

$T_{-1}=T_{1}\quad\text{and}\quad T_{n+1}=T_{n-1}\;.$ (5)

##  Appendix Appendix C Kuramoto–Sivashinsky equation

The Kuramoto–Sivashinsky equation is a nonlinear, fourth-order, one-
dimensional unsteady PDE:

$u_{t}=-(uu_{x}+u_{xx}+u_{xxxx})=-\left(\frac{1}{2}u_{x}^{2}+u_{xx}+u_{xxxx}\right)\;,$
(6)

where $u$ is the dependent chaotic variable (e.g., fluid velocity). It is
discretized similarly to the heat equation, as shown in Eq. (7), with central
differencing in space. In this case, decomposing the fourth spatial derivative
requires a five-point stencil:

$\frac{u_{i}^{m+1}-u_{i}^{m}}{\Delta
t}=-\left(\frac{(u_{i+1}^{m})^{2}-(u_{i-1}^{m})^{2}}{4\Delta
x}+\frac{u_{i+1}^{m}+u_{i-1}^{m}+2u_{i}^{m}}{\Delta x^{2}}+\right.\\\
\left.\frac{u_{i+2}^{m}-4u_{i+1}^{m}+6u_{i}^{m}-4u_{i-1}^{m}+u_{i-2}^{m}}{\Delta
x^{4}}\right)\;.$ (7)

The chaotic nature of the problem necessitates a higher-order scheme overall;
therefore, an explicit, second-order Runge–Kutta scheme, also known as the
midpoint method, is applied to the time domain.

Let the right-hand side of Eq (7) be $f(u(x,t))$, then the predictor solution
is found at $u_{i}^{m+1/2}$:

$u_{i}^{m+1/2}=u_{i}^{m}+\frac{\Delta t}{2}f(u_{i}^{m})\;.$ (8)

Then $f(u_{i}^{m+1/2})$ may be evaluated and added to $u_{i}^{m}$ to obtain
$u_{i}^{m+1}$:

$u_{i}^{m+1}=u_{i}^{m}+\Delta tf(u_{i}^{m+1/2})\;.$ (9)

The problem demonstrated here uses periodic initial and boundary conditions.
That is, the stencil at point $0$ includes points $n$ and $n-1$ and the
initial condition is periodic and continuous at the spatial boundaries.

##  Appendix Appendix D Euler equations (Sod shock tube)

The Sod shock tube problem is a one-dimensional unsteady compressible flow
problem based on the nonlinear, quasi-hyperbolic Euler equations:

$\frac{\partial Q}{\partial t}+\frac{\partial F}{\partial x}=\frac{\partial
Q}{\partial t}+J\frac{\partial Q}{\partial x}=0$ (10)

where $J$ is the Jacobian matrix,

$\displaystyle Q$ $\displaystyle=\begin{Bmatrix}\rho\\\ \rho u\\\ \rho
e\end{Bmatrix}\;,$ (11) $\displaystyle F$ $\displaystyle=\begin{Bmatrix}\rho
u\\\ \rho u^{2}+P\\\ u(\rho e+P)\end{Bmatrix}\;,$ (12)

$\rho$ is density, $e$ is internal energy, $u$ is velocity, $P$ is pressure
given by

$P=(\gamma-1)(e-\frac{\rho u^{2}}{2})\;,$ (13)

and $\gamma=1.4$ is the heat capacity ratio of air.

The initial boundary conditions, given in Eq. (14), are constant values for
the state variables on either side of a diaphragm separating two parcels of
the same fluid. The spatial boundary conditions are these values at their
respective ends of the tube:

$\begin{Bmatrix}\rho_{L}\\\ u_{L}\\\ P_{L}\end{Bmatrix}=\begin{Bmatrix}1.0\\\
0.0\\\ 1.0\end{Bmatrix}\text{ and }\begin{Bmatrix}\rho_{R}\\\ u_{R}\\\
P_{R}\end{Bmatrix}=\begin{Bmatrix}0.125\\\ 0.0\\\ 0.1\end{Bmatrix}\;.$ (14)

The equation is discretized using a second-order, finite-volume scheme with
cell-centered values. The first step in the solution is evaluating the
pressure ratio at the current timestep over a five-point stencil

$P_{r,i}=\frac{P_{i+1}-P_{i}}{P_{i}-P_{i-1}}\;\;\text{at}\;\;i-1,\;i,\;i+1$
(15)

This value is used with a minmod limiter to compute reconstructed values on
both sides, $L$ and $R$, of the current cell boundaries at $i\pm 1/2$:

$Q_{n}^{L}=\begin{cases}Q_{o}^{L}+\frac{\min(P_{r}^{L},1)}{2}*(Q_{o}^{R}-Q_{o}^{L}),&\text{if
}0<P_{r}^{L}<\infty\\\ Q_{o}^{L},&\text{otherwise.}\end{cases}$ (16)
$Q_{n}^{R}=\begin{cases}Q_{o}^{R}+\frac{\min((P_{r}^{R})^{-1},1)}{2}*(Q_{o}^{L}-Q_{o}^{R}),&\text{if
}0<\frac{1}{P_{r}^{R}}<\infty\\\ Q_{o}^{R},&\text{otherwise.}\end{cases}$ (17)

where subscript $n$ refers to the reconstructed, or new, values on the edge of
the interface and $o$ refers to the original values. For example, at $i-1/2$
the original values on the left side of the interface are at $i-1$. These
reconstructed boundaries do not represent solutions for any grid cell; they
are temporary values that interpolate the solutions.

Once we have the reconstructed values on either side of the interface, we can
calculate the flux at that cell boundary with

$\text{Flux}=\frac{1}{2}*(F(Q^{R})+F(Q^{L})+r_{sp}*(Q^{L}-Q^{R}))$ (18)

where $F(Q)$ is given by Eq. (12) and $r_{sp}$ is the spectral radius, the
largest eigenvalue of the Jacobian matrix $J$.

The spectral radius can be found with the Roe average $Q$ at the interface

$Q_{sp}=\begin{Bmatrix}\rho_{sp}\\\ u_{sp}\\\
e_{sp}\end{Bmatrix}=\begin{Bmatrix}\sqrt{\rho_{L}*\rho_{R}}\\\
\frac{\sqrt{\rho_{L}}*u_{L}+\sqrt{\rho_{R}}*u_{R}}{\sqrt{\rho_{L}}+\sqrt{\rho_{R}}}\\\
\frac{\sqrt{\rho_{L}}*e_{L}+\sqrt{\rho_{R}}*e_{R}}{\sqrt{\rho_{L}}+\sqrt{\rho_{R}}}\end{Bmatrix}$
(19)

and $P_{sp}$ with $Q_{sp}$ using Eq. (13).

The spectral radius is given by

$r_{sp}=\sqrt{\frac{\gamma*P_{sp}}{\rho_{sp}}}+|u_{sp}|$ (20)

Repeating this process at both interfaces yields all required values to solve
for a timestep

$Q_{i}^{n+1}=Q_{i}^{n}+\frac{\Delta t}{\Delta
x}(\text{Flux}_{i+1/2}^{n+1/2}-\text{Flux}_{i-1/2}^{n+1/2})$ (21)

The results presented here for the Euler equations use a second-order
Runge–Kutta scheme in time, which can be obtained with the same procedure
shown in Eqs. (9) and (8).

## Acknowledgments

This material is based upon work supported by NASA under award No. NNX15AU66A
under the technical monitoring of Drs. Eric Nielsen and Mujeeb Malik. We also
gratefully acknowledge the support of NVIDIA Corporation with the donation of
the Tesla K40c GPU used for this research.

## References

  * [1] J. P. Slotnick, A. Khodadoust, J. J. Alonso, D. L. Darmofal, W. D. Gropp, E. A. Lurie, D. J. Mavriplis, CFD vision 2030 study: A path to revolutionary computational aerosciences, NASA Technical Report, NASA/CR-2014-218178, NF1676L-18332 (Mar. 2014).
  * [2] D. A. Patterson, Latency lags bandwith, Commun. ACM 47 (10) (2004) 71–75.  doi:10.1145/1022594.1022596.
  * [3] M. Alhubail, Q. Wang, The swept rule for breaking the latency barrier in time advancing PDEs, Journal of Computational Physics 307 (2016) 110–121.  doi:10.1016/j.jcp.2015.11.026.
  * [4] M. M. Alhubail, Q. Wang, J. Williams, The swept rule for breaking the latency barrier in time advancing two-dimensional PDEs, arXiv:1602.07558 [cs.NA] (2016).
  * [5] I. Bermejo-Moreno, J. Bodart, J. Larsson, B. M. Barney, J. W. Nichols, S. Jones, Solving the compressible Navier–Stokes equations on up to 1.97 million cores and 4.1 trillion grid points, in: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis, SC ’13, ACM, New York, NY, USA, 2013, pp. 62:1–62:10.  doi:10.1145/2503210.2503265.
  * [6] K. E. Niemeyer, C. J. Sung, Recent progress and challenges in exploiting graphics processors in computational fluid dynamics, Journal of Supercomputing 67 (2) (2014) 528–564.  doi:10.1007/s11227-013-1015-7.
  * [7] NVIDIA Corporation, Whitepaper Nvidia Tesla P100, http://www.nvidia.com/object/pascal-architecture-whitepaper.html (2016).
  * [8] F. Witherden, A. Farrington, P. Vincent, PyFR: An open source framework for solving advection–diffusion type problems on streaming architectures using the flux reconstruction approach, Computer Physics Communications 185 (11) (2014) 3028–3040.  doi:10.1016/j.cpc.2014.07.011.
  * [9] M. Alhubail, Q. Wang, KSIDSwept, git commit e575d73, https://github.com/hubailmm/K-S_1D_Swept (2015).
  * [10] R. Strzodka, M. Shaheen, D. Pajak, H.-P. Seidel, Cache oblivious parallelograms in iterative stencil computations, in: Proceedings of the 24th ACM International Conference on Supercomputing, ICS ’10, ACM, New York, NY, USA, 2010, pp. 49–59.  doi:10.1145/1810085.1810096.
  * [11] T. Malas, G. Hager, H. Ltaief, H. Stengel, G. Wellein, D. Keyes, Multicore-optimized wavefront diamond blocking for optimizing stencil updates, SIAM Journal on Scientific Computing 37 (4) (2015) C439–C464.  doi:10.1137/140991133.
  * [12] K. Datta, M. Murphy, V. Volkov, S. Williams, J. Carter, L. Oliker, D. Patterson, J. Shalf, K. Yelick, Stencil computation optimization and auto-tuning on state-of-the-art multicore architectures, in: Proceedings of the 2008 ACM/IEEE Conference on Supercomputing, SC ’08, IEEE Press, Piscataway, NJ, USA, 2008, pp. 4:1–4:12.
URL http://dl.acm.org/citation.cfm?id=1413370.1413375

  * [13] M. J. Gander, 50 years of time parallel time integration, in: T. Carraro, M. Geiger, S. Körkel, R. Rannacher (Eds.), Multiple Shooting and Time Domain Decomposition Methods, Vol. 9 of Contributions in Mathematical and Computational Sciences, Springer, Cham, 2015, pp. 69–113.  doi:10.1007/978-3-319-23321-5_3.
  * [14] R. Falgout, S. Friedhoff, T. Kolev, S. MacLachlan, J. B. Schroder, Parallel time integration with multigrid, PAMM 14 (2014) 951–952.  doi:10.1002/pamm.201410456.
  * [15] M. Baboulin, S. Donfack, J. Dongarra, L. Grigori, A. Rémy, S. Tomov, A class of communication-avoiding algorithms for solving general dense linear systems on CPU/GPU parallel machines, Procedia Computer Science 9 (2012) 17–26.  doi:10.1016/j.procs.2012.04.003.
  * [16] M. Anderson, G. Ballard, J. Demmel, K. Keutzer, Communication-avoiding QR decomposition for GPUs, in: Parallel Distributed Processing Symposium (IPDPS), 2011 IEEE International, 2011, pp. 48–58.  doi:10.1109/IPDPS.2011.15.
  * [17] J. D. Owens, M. Houston, D. Luebke, S. Green, J. E. Stone, J. C. Phillips, GPU computing, Proc. IEEE 96 (5) (2008) 879–899.  doi:10.1109/JPROC.2008.917757.
  * [18] NVIDIA Corporation, CUDA C programming guide, version 8.0 (2016).
URL http://docs.nvidia.com/cuda/cuda-c-programming-guide

  * [19] A. R. Brodtkorb, T. R. Hagen, M. L. Sætra, Graphics processing unit (GPU) programming strategies and trends in GPU computing, J. Parallel Distrib. Comput. 73 (1) (2013) 4–13.  doi:10.1016/j.jpdc.2012.04.003.
  * [20] D. Storti, M. Yurtoglu, CUDA for Engineers: An introduction to High-Performance Parallel Computing, Addison-Wesley, 2015.
  * [21] F. A. Cruz, S. K. Layton, L. A. Barba, How to obtain efficient GPU kernels: An illustration using FMM & FGT algorithms, Comput. Phys. Comm. 182 (10) (2011) 2084–2098.  doi:10.1016/j.cpc.2011.05.002.
  * [22] D. J. Magee, K. E. Niemeyer, Data, plotting scripts, and figures for “Accelerating solutions of PDEs with GPU-based swept time-space decomposition” (Nov. 2017).  doi:10.6084/m9.figshare.4968050.v4.
  * [23] Q. Wang, Decomposition of stencil update formula into atomic stages, arXiv:1606.00721 [math.NA] (2017).
  * [24] M. Harris, CUDA pro tip: Do the Kepler shuffle, https://devblogs.nvidia.com/parallelforall/cuda-pro-tip-kepler-shuffle/, accessed: 3 June 2016 (Feb. 2014).
  * [25] D. J. Magee, K. E. Niemeyer, Niemeyer-Research-Group/1DSweptCUDA v2 (May 2017).  doi:10.5281/zenodo.570984.
  * [26] S. Xiao, W. C. Feng, Inter-block GPU communication via fast barrier synchronization, in: 2010 IEEE International Symposium on Parallel Distributed Processing (IPDPS), 2010, pp. 1–12.  doi:10.1109/IPDPS.2010.5470477.


#  Monte Carlo Study of Relaxor Systems:
A Minimum Model of Pb(In1/2Nb1/2)O3

Yusuke Tomita E-mail address<EMAIL_ADDRESS>Takeo Kato and Kazuma
Hirota1 Institute for Solid State PhysicsInstitute for Solid State Physics
University of Tokyo University of Tokyo Kashiwa 277-8581 Kashiwa 277-8581
Japan
1Department of Earth and Space Science Japan
1Department of Earth and Space Science Faculty of Science Faculty of Science
Osaka University Osaka University Toyonaka 560-0043 Toyonaka 560-0043 Japan
Japan

###### Abstract

We examine a simple model of Pb(In1/2Nb1/2)O3 (PIN), which includes both long-
range dipole-dipole interaction and random local anisotropy. An improved
algorithm optimized for long-range interaction has been applied to efficient
large-scale Monte Carlo simulation. We demonstrate that the phase diagram of
PIN is qualitatively reproduced by this minimum model. Some characteristic
features of relaxors such as nanoscale domain formation, slow dynamics, and
dispersive dielectric responses are also examined.

relaxors, ferroelectric phase transition, antiferroelectric phase transition,
Monte Carlo algorithm, domain formation, slow dynamics, dispersive dielectric
responses, Pb(In1/2Nb1/2)O3

Relaxors, whose first discovery dates back to over half a century ago [1],
have newly attracted much interest over the past few decades owing to their
colossal dielectric and piezoelectric responses that are appealing for
industrial applications. Despite intensive research, however, the physical
origin of the unusual properties of relaxors is not fully understood yet. [2]
This is because nanoscale intrinsic randomness in relaxor systems has to be
addressed appropriately.

As a simple example to understand the physics of relaxors, let us choose the
lead-based relaxor Pb(In1/2Nb1/2)O3 (PIN), which is the target of the present
theoretical work. The intrinsic inhomogeneities of relaxor systems originate
from the random configuration of In3+ and Nb5+ at the B-site in the perovskite
structure, while the dielectric property of such systems is governed mainly by
the replacement of Pb and O ions. This system has the advantage that the
strength of randomness can be controlled by adjusting annealing temperature.
[3] For the lowest annealing temperature, an alternate order of In and Nb
atoms stabilizes the fourfold antiferroelectric (AFE) phase. For the higher
annealing temperature, an increase in B-site randomness decreases AFE
transition temperature. Under sufficiently strong randomness in the B-site,
the ferroelectric phase accompanies relaxor properties. This behavior of PIN
suggests that the potential FE instability exists behind the AFE phase, and
that its emergence by the selective suppression of the AFE phase due to B-site
randomness is relevant to FE relaxors. Recent experimental progress in the
X-ray diffraction analysis facilitated by strong-intensity photon sources has
also enabled us to detect the potential competition between the FE and AFE
phases of PIN from the viewpoint of lattice dynamics. [4, 5]

There are several theoretical methods for elucidating the emergence of relaxor
properties. They utilize the existing effective theory dealing with a strong
randomness such as the extension of the Ginzburg-Landau-Devonshire theory, [6,
7] the spherical random-bond-random-field model, [8] and the dynamical-related
model. [9, 10] Although these approaches have provided useful description of
relaxor properties, it is unclear how their model parameters should be derived
microscopically. At present, first-principles calculation is playing a crucial
role in the microscopic description of ferroelectrics and ferroelectric
relaxors. In the pioneering work by Cohen and coworker [11, 12], the electron-
structure calculation of BaTiO3 and PbTiO3 indicated the importance of the
covalent nature among composite atoms through Ti-O hybridization. After the
success of their work, microscopic first-principles approaches have been
adopted in a number of issues. Further theoretical progress has been achieved
by a hybrid method combining the microscopic derivation of the effective
Hamiltonian and its numerical simulation. Zhong et al. have derived the
effective Hamiltonian of BaTiO3 by extracting the adiabatic potential of
important phonon modes from first-principles calculation. [13, 14] Monte Carlo
simulation for this effective Hamiltonian has enabled the successful
reproduction of the sequential finite-temperature phase transitions of BaTiO3
in good agreement with experimental results. Recently, the hybrid method has
also been applied to the study of relaxor systems. [15, 16]

Thus, the utilization of first-principles calculation is a promising method
for understanding relaxors. However, there still remain several difficulties
to overcome for the actual application of this method to the study of
relaxors. One major difficulty is in the use of the numerical solver of the
effective Hamiltonian for phonons. Even though the effective Hamiltonian
derived for relaxors is in a simple form, its numerical simulation may become
extremely difficult because it is inevitable to encounter slow dynamics
inherent to random systems. The Monte Carlo simulation of the effective
Hamiltonian needs a large number of Monte Carlo steps owing to its long
correlation time. Therefore, the development of appropriate numerical solvers
for relaxor systems seems indispensable. This situation may remind us of Monte
Carlo studies of spin-glass systems. Although the Hamiltonian of spin-glass
systems has a simple form, much effort has to be exerted for sufficient Monte
Carlo sampling to understand their peculiar slow dynamics. In the research
field of relaxors, however, the importance of developing a numerical solver
has not been addressed so far.

In this paper, we propose a simple model of ferroelectric relaxors to explain
the phase diagram of PIN. Our model includes both long-range dipole-dipole
interaction and local randomness. We apply the new effective algorithm
optimized for long-range interaction proposed by Fukui and Todo. [17] This new
algorithm enables us to simulate long-range interaction systems with the cost
of O($N$) with respect to the system size $N$, while the conventional
simulation for the same system takes the cost of O($N^{2}$). In addition to
the substantial reduction in computational time by employing this algorithm,
we adopt the exchange Monte Carlo method to attenuate the slow dynamics of
random systems. [18] We show that our simple model may qualitatively reproduce
the phase diagram of PIN. Moreover, we demonstrate that some characteristic
features of relaxors such as domain formation and dispersive dielectric
response are reproduced reasonably. Our results indicate the power of the
sophisticated Monte Carlo method as well as the possibility that relaxor
systems may be represented by a simple model.

We consider the model Hamiltonian for PIN on a 2D square lattice as

$\displaystyle{\cal H}$ $\displaystyle=$
$\displaystyle\sum_{i<j}\left[\frac{\mbox{\boldmath$S$}_{i}\cdot\mbox{\boldmath$S$}_{j}}{r^{3}_{ij}}-3\frac{(\mbox{\boldmath$r$}_{ij}\cdot\mbox{\boldmath$S$}_{i})(\mbox{\boldmath$r$}_{ij}\cdot\mbox{\boldmath$S$}_{j})}{r^{5}_{ij}}\right]$
(1) $\displaystyle-$
$\displaystyle\sum_{i}(\mbox{\boldmath$D$}_{i}\cdot\mbox{\boldmath$S$}_{i})^{2},$

where $\mbox{\boldmath$S$}_{i}$ is a unit vector in the $xy$-plane
representing the dipole moment on the $i$th unit cell induced by the off-
center replacement of the Pb atom. The first term of the Hamiltonian is the
dipole-dipole interaction dependent on the relative position
$\mbox{\boldmath$r$}_{ij}=\mbox{\boldmath$r$}_{j}-\mbox{\boldmath$r$}_{i}$
between the sites $i$ and $j$, while the second term describes local
anisotropy whose direction and strength are denoted as
$\mbox{\boldmath$D$}_{i}$. In order to reproduce the phase diagram of PIN, we
design our model such that the FE phase is stabilized by dipole-dipole
interaction, while the AFE phase is stabilized in the presence of an
alternative change in $\mbox{\boldmath$D$}_{i}$. Supposing that B-site
randomness affects only the local energy change through the anisotropy
parameter $\mbox{\boldmath$D$}_{i}$, we can expect that all the features of
PIN, i.e., the AFE transition in the ordered PIN, its suppression by B-site
randomness, and the appearance of the FE domain in the sufficiently disordered
PIN are reproduced.

The detailed setting of our model is given as follows. In a naive 2D square
lattice, the dipole-dipole interaction does not lead to ferroelectric
instability because the columnar antiferroelectric state is more favored in
orthogonal lattices. [19] Therefore, we modify our model slightly: We divide
the unit cells into two bipartite sublattices named P and Q, and shift only
the Q-sites in the $z$-direction by a unit length. This rearrangement of unit
cells ensures ferroelectric instability in the absence of local anisotropy.
The antiferroelectric instability is driven by an alternative arrangement of
two types of anisotropy, namely, $\mbox{\boldmath$D$}_{1}$ and
$\mbox{\boldmath$D$}_{2}$ in the P- and Q-sites, where
$\mbox{\boldmath$D$}_{1}=(D/\sqrt{2},-D/\sqrt{2})$ and
$\mbox{\boldmath$D$}_{2}=(D/\sqrt{2},D/\sqrt{2})$. The randomness of $D$ is
controlled by the probability $p$, where $\mbox{\boldmath$D$}_{1}$
($\mbox{\boldmath$D$}_{2}$) is attached to the P-(Q-)site with a probability
$(1-p)^{2}$, oppositely to $p^{2}$, and is turned off ($D$=0) with $2p(1-p)$.
The two limiting values $p=0$ and $p=1/2$ correspond to the ordered and
completely disordered PINs, respectively.

The model Hamiltonian is examined by Monte Carlo simulations under a periodic
boundary condition to minimize the strong surface effect due to long-range
interaction. Long-range interactions are summed up by the Ewald summation
technique. [20] We employ the ${\rm O}(N)$ method based on Walker’s algorithm
for efficient update [17, 21] as well as the temperature-exchange algorithm.
[18] One exchange trial between replicas was made for each of the $10$ MC
steps. The system size is taken up to $N=32\times 32$, for which $3\times
10^{6}$ MC steps for thermalization and $2\times 10^{5}$ MC steps for
measurement were needed in the severest case of $p=1/2$. The sample average is
taken over $10$ different random configurations of $\mbox{\boldmath$D$}_{i}$
for each $p$. Throughout this paper, the strength of the anisotropy is fixed
as $|\mbox{\boldmath$D$}_{i}|=1$.

Figure 1:  (Color online) (a) A plot of four-fold antiferroelectric(AFE) order
parameters in the absence of randomness ($p=0$) and (b) a plot of
ferroelectric(FE) order parameters at $p=0.5$ are shown for different system
sizes $N=L\times L$. (c) Phase diagram obtained by the minimum model. Phase
boundaries in the figure are visual guides.

Figure 1(a) shows the temperature dependence of the average squared four-fold
staggered polarization $m_{\rm AF}$ for the ordered PIN ($p=0$). The pattern
of the AFE ordering, which is shown in Fig. 1(c), agrees with experimental
observations. As the system size $N=L\times L$ increases, a sharp increase in
the order parameter appears below a critical temperature, indicating the phase
transition. The transition point is determined by the crossing point of the
Binder parameters, $\langle m_{\rm AF}^{4}\rangle/\langle m_{\rm
AF}^{2}\rangle^{2}$, of $L=16$ and $32$. By a similar analysis of nonzero
values of $p$, the phase boundary of AFE in the $p$-$T$ plane is determined,
as shown by full circles in Fig. 1(c). For small values of $p$, the transition
temperature of AFE is suppressed by the B-site randomness. For sufficiently
large values of $p$, the FE domain develops at low temperatures instead of the
AFE phase. Figure 1(b) shows the average squared uniform polarization as a
function of temperature in the completely disordered case ($p=1/2$). An abrupt
increase in squared polarization below a threshold temperature for the
$N=32\times 32$ system indicates a rapid development of the FE domain. In our
simulation, however, no long-range FE ordering could be detected by finite-
size scaling because of the very slow relaxation of the Monte Carlo sampling
and the complex size dependence in the presence of strong randomness. Here, we
roughly estimate the threshold temperature at which the local FE instability
rises up by Binder-parameter analysis between $L=16$ and $32$, and plot it
using the empty circle in Fig. 1(c). At approximately $p=0.4$, the competition
between AFE and FE makes the Monte Carlo calculation severe, and the phase
boundary could not be determined. Rough interpolations for AFE and FE are
shown by the solid curves in Fig. 1(c) only as visual guides. The obtained
phase diagram is in good agreement with experimental results of PIN.

Figure 2:  (Color online) Log-log plot of the normalized time autocorrelation
functions. While the autocorrelation functions decay smoothly for the no-
anisotropy system at $T=0.2$ and the random system at $T=0.3$, the
autocorrelation function shows a slow decay for the random system at $T=0.2$.

For large values of $p$, some characteristic features of relaxor systems
appear below a crossover temperature, which locates higher than the threshold
temperature for FE domain formation, as indicated roughly by the dashed line
in Fig. 1(c). The emergence of relaxor properties is indicated significantly
by the slow convergence of the Monte Carlo average. Figure 2 shows the
normalized autocorrelation of the energy $E$ for $N=32\times 32$ defined by

$\phi(t)=\frac{\langle E(0)E(t)\rangle-\langle E\rangle^{2}}{\langle
E^{2}\rangle-\langle E\rangle^{2}}$ (2)

as a function of the Monte Carlo step $t$ in the completely random case
($p=1/2$). Here, for the purpose of examining the dynamics of the system, we
turned off the exchange process between replicas. We find that the decay of
the autocorrelation function is extremely slow at $T=0.2$ in clear contrast
with that at a higher temperature $T=0.3$. This slow energy relaxation
indicates a glasslike behavior for relaxors at low temperatures. In order to
confirm that this slow relaxation comes from randomness of the B-site, we also
calculated the autocorrelation at $T=0.2$ in the absence of local anisotropy
($|\mbox{\boldmath$D$}_{i}|=0$), for which normal FE ordering is expected. As
seen in Fig. 2, the autocorrelation function follows a simple decay with no
anomalous slow relaxation. Thus, the present model is expected to exhibit
relaxor properties under sufficiently strong randomness below a crossover
temperature.

Figure 3:  (Color online) Snapshots of dipole configuration for $p=0.5$ and
$T=0.1$. Each color denotes the angle of dipole moment. Reference color bars
are shown on the right. The bottom color denotes $-\pi$, and the top color
denotes $\pi$. The last graph shows the result for the same parameters but for
a model with a dipole interaction interrupted up to the next-nearest sites.

In order to visualize the glassy state realized in the relaxor phase, we show
in Fig. 3 the spatial pattern of dipole directions by taking snapshots as a
function of the MC steps for $p=0.5$ and $N=32\times 32$. We can see that
mosaic-like FE domains are formed after a finite relaxation time. The
boundaries (domain-wall regions) between the neighboring FE domains are rather
clear. One may consider that the FE domain structure is determined by a
partially ordered configuration of anisotropy, the so-called chemical
nanoregion (CNR). The present FE domain is, however, irrelevant to CNR,
because partial ordering is absent for a complete random choice of anisotropy
at each site in our model.

As has already been mentioned, the present system possesses FE instability in
the absence of local anisotropy (${\boldmath D}={\boldmath 0}$). The random
configuration of anisotropy similarly causes a local FE instability due to
effective cancellation of anisotropy, but at the same time prevents the growth
of an FE domain into a uniform FE order, as seen in Fig. 3. A small but finite
change of the snapshot in the long-time region of the Monte Carlo simulation
indicates that a very slow fluctuation of the FE domain governs the relaxor
properties in the present model. Here, we should note that the evolution of
the dipole configuration along the Monte Carlo step is not directly related to
actual real-time dynamics. The snapshots, however, provide intuitive and
fruitful insights for understanding physical processes in relaxor systems. In
our model, the long-range nature of the dipole interaction is essential for
the formation of the FE domain structure. To see this, we show a snapshot of
Monte Carlo simulation for a model with dipole interaction interrupted up to
the next-nearest site in the last graph of Fig. 3. We find that the domain
structure disappears in this modified model.

Figure 4: (Color online) Dielectric constants evaluated as functions of
temperature under the ac electric field for several frequencies in the (a)
ordered case ($p=0$) and (b) completely disordered case ($p=1/2$).

The existence of nanoscale frustrated FE domains indicates the high degeneracy
of glassy states. This degeneracy is expected to be responsible for large
dielectric responses since a significant change in the polarization of such FE
domains is possible by a small external electric field. To see this, we
measure the Monte Carlo evolution of the total polarization under a small ‘ac’
electric field, which varies periodically along the Monte Carlo step. We then
calculate the ‘ac’ dielectric constant in the linear response regime. We
should note again that the ‘ac’ dielectric constant thus obtained is not equal
to the real-time ac dielectric constant. Furthermore, we should keep in mind
that this is a nonequilibrium response since the relaxation to equilibrium
states is never realized in the glassy phase in the present simulation.
Nevertheless, the ‘ac’ dielectric response calculated here is expected to
mimic the actual ac response in a qualitative level. In Fig. 4, we show the
‘ac’ dielectric constant as a function of temperature for three frequencies in
the ordered ($p=0.0$) and completely disordered ($p=1/2$) cases. In both
cases, dielectric constant shows a peak at approximately the transition
(crossover) temperature. The marked difference is found in the low-temperature
phase. In the ordered case, dielectric constant sharply drops in the low-
temperature AFE phase, and becomes almost independent of frequency. On the
other hand, in the disordered case, it decreases gradually as temperature
decreases, and a strong frequency dependence remains. This strong dispersion
of dielectric constant at low temperatures can be explained as follows. Under
a high-frequency electric field, each dipole inside FE domains may respond to
an external field. Under a low-frequency electric field, however, those with
frustration that construct large domains start to respond. Therefore, in a
minimum model, a dipolar glass with ferroelectric ordering is realized. This
ordering causes a broadened dielectric constant and a strong dependence on
frequency in the relaxor phase.

In summary, we examined a simple theoretical model of PIN composed of dipolar
interaction and local random anisotropy. We demonstrated that efficient Monte
Carlo simulations equipped with an improved algorithm optimized for long-range
interaction may access several characteristic features of relaxors. The phase
diagram of PIN was qualitatively reproduced by appropriate inclusion of the
intrinsic competition between the AFE and FE phases. By the examination of the
Monte Carlo evolution of the dipole configuration, we demonstrated some of the
glassy behaviors inherent to relaxor systems such as the FE domain formation,
extremely slow dynamics, and strong frequency dependence of dielectric
responses.

We end this paper by mentioning the future outlook of theoretical approaches
to studying relaxors. The model we treated here includes a few artificial
assumptions. They, however, will be removed by replacing our model with the
effective Hamiltonian derived by first-principles calculation. We stress that
the smart Monte Carlo algorithm is applicable not only to the rotator model
but also to the continuous-variable model. The hybridization of the first-
principles calculation and statistical approach based on the Monte Carlo
simulation will be an effective means of elucidating the microscopic origin of
relaxors.

The computation in the present work is executed on computers at the
Supercomputer Center, Institute for Solid State Physics, University of Tokyo.
The present work is financially supported by a MEXT Grant-in-Aid for
Scientific Research (B) (19340109), by a MEXT Grant-in-Aid for Scientific
Research on Priority Areas “Novel States of Matter Induced by Frustration”
(19052002), and by the Next Generation Supercomputing Project, Nanoscience
Program, MEXT, Japan

## References

  * [1] G. A. Smolenskii and A. I. Agronovskaya: Sov. Phys. Tech. Phys. A3 (1958) 1380.
  * [2] For recent reviews, A. A. Bokov and Z.-G. Ye: J. Mater. Sci. 41 (2006) 31.
  * [3] A. A. Bokov, I. P. Raevskii, O. I. Prokopalo, E. G. Fesenko, and V. G. Smotrakov: Ferroelectrics 54 (1984) 241.
  * [4] K. Ohwada, K. Hirota, H. Terauchi, T. Fukuda, S. Tsutui, A. Q. R. Baron, J. Mizuki, H. Ohwa, and N. Yasuda: Phys. Rev. B 77 (2008) 094136.
  * [5] K. Ohwada and Y. Tomita: J. Phys. Soc. Jpn. 79 (2010) 011012.
  * [6] A. F. Devonshire: Adv. Phys. 3 (1954) 85.
  * [7] J.-M. Liu, S. T. Lau, H. L. W. Chan, and C. L. Choy: J. Mater. Sci. 41 (2006) 163.
  * [8] R. Blinc, J. Dolinšek, A. Gregorovič, B. Zalar, C. Filipič, Z. Kutnjak, A. Levstik, and R. Pirc: Phys. Rev. Lett. 83 (1999) 424.
  * [9] A. J. Bell: J. Phys.: Condens. Matter 5 (1993) 8773.
  * [10] B. E. Vugmeister and H. Rabitz: Phys. Rev. B 57 (1998) 7581.
  * [11] R. E. Cohen: Nature 358 (1992) 136.
  * [12] R. E. Cohen and H. Krakauer: Ferroelectrics 136 (1992) 65.
  * [13] W. Zhong, D. Vanderbilt, and K. M. Rabe: Phys. Rev. Lett. 73 (1994) 1861.
  * [14] W. Zhong, D. Vanderbilt, and K. M. Rabe: Phys. Rev. B 52 (1995) 6301.
  * [15] S. Tinte, B. P. Burton, E. Cockayne, and U. V. Waghmare: Phys. Rev. Lett. 97 (2006) 137601.
  * [16] B. P. Burton, E. Cockayne, S. Tinte, and U. V. Waghmare: Phys. Rev. B 77 (2008) 144114.
  * [17] K. Fukui and S. Todo: J. Comp. Phys. 228 (2009) 2629.
  * [18] K. Hukushima and K. Nemoto: J. Phys. Soc. Jpn. 65 (1996) 1604.
  * [19] Y. Tomita: J. Phys. Soc. Jpn. 78 (2009) 114004.
  * [20] J.-J. Weis: J. Phys.: Condens. Matter 15 (2003) S1471.
  * [21] A. J. Walker: ACM Trans. Math. Software 3 (1977) 253.