Analysis of heterogeneous computing approaches to simulating heat transfer in heterogeneous material

TL;DR

This paper presents a parallelizable, assembly-free finite element method for simulating heat transfer in heterogeneous materials, optimized for GPU computing, and demonstrates its efficiency on large-scale, real-world problems.

Contribution

The paper introduces a novel assembly-free FEM approach adapted for GPU architectures, enabling efficient simulation of heat transfer in complex heterogeneous materials.

Findings

GPU implementation outperforms CPU sparse matrix methods

Dual GPU setup accelerates simulations significantly

Method effectively handles large-scale, real-world inverse problems

Abstract

The simulation of heat flow through heterogeneous material is important for the design of structural and electronic components. Classical analytical solutions to the heat equation PDE are not known for many such domains, even those having simple geometries. The finite element method can provide approximations to a weak form continuum solution, with increasing accuracy as the number of degrees of freedom in the model increases. This comes at a cost of increased memory usage and computation time; even when taking advantage of sparse matrix techniques for the finite element system matrix. We summarize recent approaches in solving problems in structural mechanics and steady state heat conduction which do not require the explicit assembly of any system matrices, and adapt them to a method for solving the time-depended flow of heat. These approaches are highly parallelizable, and can be performed on graphical processing units (GPUs). Furthermore, they lend themselves to the simulation of heterogeneous material, with a minimum of added complexity. We present the mathematical framework of assembly-free FEM approaches, through which we summarize the benefits of GPU computation. We discuss our implementation using the OpenCL computing framework, and show how it is further adapted for use on multiple GPUs. We compare the performance of single and dual GPUs implementations of our method with previous GPU computing strategies from the literature and a CPU sparse matrix approach. The utility of the novel method is demonstrated through the solution of a real-world coefficient inverse problem that requires thousands of transient heat flow simulations, each of which involves solving a 1 million degree of freedom linear system over hundreds of time steps.

Tables

Table 1. Table 1: Ratio of the rate of increase of wall clock time with increasing degrees of freedom between the sparse matrix CPU method and the linear asymptote of each GPU method.

Method	1	2	3: single GPU	3: dual GPU
Speedup Factor	6.5	1.2	8.4	16.8

Table 2. Table 2: Full performance comparison between the sparse matrix CPU method with incomplete Cholesky preconditioning and our best performing implementation.

Degrees of Freedom	10.5$\times {10}^3$	78.1$\times {10}^3$	257$\times {10}^3$	600$\times {10}^3$	1.16$\times {10}^6$	2.00$\times {10}^6$
Sparse matrix CPU: Total time (s)	0.2	1.6	6.2	16	37	67
Total Iterations	39	47	54	56	62	67
Time per Iteration (s)	5.6$\times {10}^{-3}$	34$\times {10}^{-3}$	110$\times {10}^{-3}$	280$\times {10}^{-3}$	600$\times {10}^{-3}$	1000$\times {10}^{-3}$
FG DbD: Total time (s)	0.48	1.7	2.0	5.4	13	26
Total Iterations	287	344	567	780	1047	1278
Time per Iteration (s)	1.6 $\times {10}^{-3}$	5$\times {10}^{-3}$	3.6$\times {10}^{-3}$	6.9$\times {10}^{-3}$	12$\times {10}^{-3}$	20$\times {10}^{-3}$

Table 3. Table 3: Memory buffer types and initialization values. If a number is given for initialization, the specified number of bytes is allocated.

Name	Flag	Initialization
M_FG_buf	HOST_TO_DEVICE_USE	M_FG
K_FG_buf	HOST_TO_DEVICE_USE	K_FG
DoFMapLocal_buf	HOST_TO_DEVICE_COPY	DoFMapLocal
coordinateMapLocal_buf	HOST_TO_DEVICE_COPY	coordinateMapLocal
C_buf	HOST_TO_DEVICE_COPY	C
vert_scale_buf	HOST_TO_DEVICE_COPY	vert_scale
corr_bounds_buf	HOST_TO_DEVICE_COPY	corr_bounds
mat_coefs_buf	HOST_TO_DEVICE_COPY	mat_coefs
P_buf	DEVICE_READ_WRITE	x.nbytes
u0_buf	HOST_TO_DEVICE_COPY	u0
Fdt_buf	HOST_TO_DEVICE_COPY	Fdt
x_buf	HOST_READ_WRITE	x
b_buf	DEVICE_READ_WRITE	x.nbytes
r_buf	DEVICE_READ_WRITE	x.nbytes
d_buf	DEVICE_READ_WRITE	x.nbytes
q_buf	DEVICE_READ_WRITE	x.nbytes
s_buf	DEVICE_READ_WRITE	x.nbytes
delta_buf	PINNED	delta
alpha_buf	DEVICE_READ_WRITE	delta.nbytes
neg_alpha_buf	DEVICE_READ_WRITE	delta.nbytes
delta_new_buf	DEVICE_READ_WRITE	delta.nbytes
beta_buf	DEVICE_READ_WRITE	delta.nbytes
r1_buf	DEVICE_READ_WRITE	r1_size$\times$4
r2_buf	DEVICE_READ_WRITE	r2_size$\times$4
VVM_loc_buf	LocalMemory	max_wg_size$\times$4
Ax_split_buf	DEVICE_READ_WRITE	x.nbytes$\times$4
P_split_buf	DEVICE_READ_WRITE	x.nbytes$\times$4
x_local_buf	LocalMemory	32$\times$4$\times$4

Argument	In/Out	Memory space	Data type	Size
*x	input	global	float64	nVertices
*y	input	global	float64	nVertices
nVertices	input		uint32	1
*VVM_loc		local	float64	max_workgroup_size
*r	output	global	float64	r1_size

Argument	In/Out	Memory space	Data type	Size
*rPrevious	input	global	float64	r(k-1)_size $\times$ max_workgroup_size
nVertices	input		uint32	1
*VVM_loc		local	float64	max_workgroup_size
*r	output	global	float64	rk_size

Argument	In/Out	Memory space	Data type	Size
*rPrevious	input	global	float64	r(k-1)_size $\times$ max_workgroup_size
*delta	input	global	float64	1
nVertices	input		uint32	1
*b		local	float64	max_workgroup_size
*alpha	output	global	float64	1
*neg_alpha	output	global	float64	1

Argument	In/Out	Memory space	Data type	Size
*x	input	global	float64	nVertices
*y	input	global	float64	nVertices
*a	input	global	float64	1
*x_plus_ay	output	global	float64	nVertices

Argument	In/Out	Memory space	Data type	Size
*P	input	global	float64	nVertices
*x	input	global	float64	nVertices
*Pinvx	output	global	float64	nVertices

Argument	In/Out	Memory space	Data type	Size
*delta_new	input	global	float64	1
*delta_old	both	global	float64	1
*beta	output	global	float64	1

Argument	In/Out	Memory space	Data type	Size
*u0	both	global	float64	nVertices
*u_new	input	global	float64	nVertices

Argument	In/Out	Memory space	Data type	Size
*M_FG	input	constant	float64	4
*K_FG	input	constant	float64	4
*x	input	global	float64	nVertices
*DoFMapLocal	input	constnat	uint32	12
*coordinateMapLocal	input	constant	float64	12
*C	input	constant	uint32	3
*vert_scale	input	constant	float64	6
*corr_bounds	input	constant	float64
*mat_coefs	input	constant	float64	4
theta	input		float64	1
dt	input		float64	1
*x_local		local	float64	4$\times$32
*Ax_split_local		local	float64	4$\times$32$\times$12
*Ax_split	output	global	float64	nVertices$\times$4

Argument	In/Out	Memory space	Data type	Size
*Ax_split	input	global	float64	nVertices$\times$4
*Ax	output	global	float64	nVertices

Argument	In/Out	Memory space	Data type	Size
*Ax_split	input	global	float64	nVertices$\times$4
*b	input	global	float64	nVertices
c	input		float64	1
*cAx_plus_b	output	global	float64	nVertices

Equations

$$\{\begin{array}{ll}\rho C\frac{\partial T(\vec{x},t)}{\partial t}=\nabla \cdot (k\nabla T(\vec{x},t))& \text{mbox}indomain\Omega \times (0,t_f),\\ k\frac{\partial T(\vec{x},t)}{\partial \vec{n}}=f(\vec{x})& \text{mbox}onsides,\\ T(\vec{x},0)=T_{\text{ambient}},& \end{array}$$

$$a(T^{(i)}(\vec{x}),\varphi (\vec{x}))$$

$$L(\varphi (\vec{x}))$$

$$\mathbf{M}={\left[{\int }_{\Omega }\rho C{\varphi }_m{\hat{\varphi }}_{n}d\vec{x}\right]}_{m,n\in \{1,\dots ,N\}}\text{ and }\mathbf{K}={\left[{\int }_{\Omega }k\nabla {\varphi }_m\cdot \nabla {\hat{\varphi }}_{n}d\vec{x}\right]}_{m,n\in \{1,\dots ,N\}}.$$

$$\left[\mathbf{M}+\theta \Delta t\mathbf{K}\right]{\vec{U}}^{(i)}=\left[\mathbf{M}-(1-\theta )\Delta t\mathbf{K}\right]{\vec{U}}^{(i-1)}+\Delta t\vec{F}$$

$${\mathbf{M}}_e={\left[{\int }_{{\Omega }_e}{\varphi }_m{\hat{\varphi }}_{n}d\vec{x}\right]}_{m,n\in \{1,\dots ,D\}}.$$

$$\text{vbox}\mathsf{A}{\text{displaylimits}}_{e\in \mathcal{E}}(\rho C{)}_e{\mathbf{M}}_e=\mathbf{M}.$$

$$\text{vbox}\mathsf{A}{\text{displaylimits}}_{e\in \mathcal{E}}(\rho C{)}_e{\mathbf{M}}_e{\vec{x}}_e=\text{vbox}\mathsf{A}{\text{displaylimits}}_{n\in \mathcal{N}}\left(\sum _{e\in {\mathcal{E}}^{(n)}}(\rho C{)}_e{\mathbf{M}}_e^n{\vec{x}}_e\right)=\vec{y},$$

$$\mathbf{A}{\vec{U}}^i=\text{vbox}\mathsf{A}{\text{displaylimits}}_{e\in \mathcal{E}}{\mathbf{A}}_e{\vec{U}}_e^i=\text{vbox}\mathsf{A}{\text{displaylimits}}_{n\in \mathcal{N}}\left(\sum _{e\in {\mathcal{E}}^{(n)}}{\mathbf{A}}_e^n{\vec{U}}_e^i\right)=\vec{b}.$$

$${\vec{y}}_i=\left[\sum _{\begin{array}{c}e\in {\mathcal{E}}^{(i)}\\ i={\mathcal{N}}^{(e)}(j)\end{array}}\left[{\mathbf{A}}_e^j{\vec{x}}_e\right]\right](i\in \mathcal{N}),$$

$${\vec{y}}_i=\left[\sum _{\begin{array}{c}e\in {\mathcal{E}}^{(i)}\\ i={\mathcal{N}}^{(e)}(j)\end{array}}{\mathbf{A}}_e^j{\vec{x}}_e\right](i\in \mathcal{N}).$$

$$\rho C=\{\begin{array}{ll}3.724\text{e}6{\text{ g/mm C s}}^2& x_3\le 5\\ 1.65\text{e}6{\text{ g/mmCs}}^2& x_3>5\end{array},k=\{\begin{array}{ll}4.9\text{e}8{\text{ mm}}^2{\text{C s}}^2& x_3\le 5\\ 4\text{e}6{\text{ mm}}^2{\text{C s}}^2& x_3>5\end{array},$$

$${\theta }_{i+1}=\{\begin{array}{ll}\hat{\theta }& \text{ with probability }\min (1,p(\mathcal{D}|\hat{\theta })/p(\mathcal{D}|{\theta }_i))\\ {\theta }_i& \text{ otherwise.}\end{array}$$

$$\text{mbox}r1\_size=\lceil \frac{\text{mbox}nVertices}{\text{mbox}max\_workgroup\_size}\rceil .$$

$$\text{mbox}rk\_size=\lceil \frac{\text{mbox}r(k-1)\_size}{\text{mbox}max\_workgroup\_size}\rceil .$$

$$\text{mbox}MVM\_global\_size=186\lceil \frac{6(C[0]+1)(C[1]+1)C[2]}{180}\rceil .$$

Full text

1 Introduction

The finite element method (FEM) provides useful approximations to weak solutions of partial differential equations (PDEs) over a much broader family of domains than those for which analytical solutions are known [19]. The method involves transforming a continuum solution into the solution of a sparse system of linear equations, which are well understood and can be quickly solved or approximated. This transformation encodes the domain geometry as well as properties of the domain that can vary over space and time, such as material properties for heat conduction problems. Fast solvers for this type of problem are critical for providing simulations over domains with increasingly fine discretization. This is particularly salient in the case that many successive simulations are desired with changes to the parameterization of material coefficients, for example: as in the solution of a coefficient inverse problem [6].

The present study is motivated by such inverse problems. Furthermore, increasingly efficient methods for the simulation of heat transfer through heterogeneous material are useful for the design of structural and electrical systems [14]. In this area of research, analytical solutions have been found in some restricted cases [32], and only some characteristics of the transient behavior of more general problems are known [13, 1, 29]. The purely mathematical analysis relies on isotropy or bulk approximation over the domain. The presence of more than one set of heat conduction parameters over a heterogenous material makes the problem too complex to solve, or results in solutions that are too complicated to be useful [13].

1.1 Scope and Organization

This paper is divided into six sections. In Section 2, some history of matrix-free methods for FEM is discussed; leading to recent heterogeneous computing approaches. In Section 3, the necessary mathematical preliminaries for our methods are summarized. The assembly operator is introduced, along with the flexibility that it offers for matrix-free methods. This flexibility is explored in Section 4 with three interpretations of the assembly operator and the resulting implementations for the simulation of heat transfer through heterogeneous media. In Section 5, numerical experiments and their results are described, to compare the speed of the implementations over a range of problem sizes. A sparse serial method is also included in the comparison. A coefficient inverse problem arising from a real-world need for nondestructive corrosion detection is introduced and analyzed in Section 5.4. Finally, Section 6 summarizes the goals and future impact of this research.

2 Background and Motivation

The methods that are developed in this work are derived from an element-by-element (EbE) decomposition of the finite element method. This viewpoint was introduced by Hughes et al. in 1983, for heat conduction calculations that were large for the time [12]. The motivation at that time was in avoiding exceedance of available computer memory for the storage of the large FEM system matrix, rather than considerations of speed. The authors showed stability for the algorithm, and followed up to show the same properties of the EbE method for structural and solid mechanics problems [11]. Carey et al. exposed the potential for the EbE method to be parallelized, with a demonstration of a 2D convection-diffusion simulation in 1988 [3]. Following this, research interest in the EbE framework had little room for advancement until advances in computer hardware were made.

Heterogeneous computing is the practice of using dissimilar coprocessors in the solution of a numerical problem; typically a computer’s central processing unit (CPU) and one or more graphics processing units (GPUs). This provides the user with access to fast serial processing power, as well as large-scale parallelism for certain computations that are properly amenable. Kiss et al. first utilized GPUs for the EbE method, using NVIDIA’s proprietary CUDA platform [15]. The authors explored considerations that are particular to GPU computing, such as the preference for repeated computation rather than loading data from memory, and segmenting the domain with graph-coloring methods to avoid race conditions in parallelism. There have since been other efforts to parallelize PDE solvers for GPUs that are leveraged on symmetries of a particular problem [25]. A noteable recent package is the TeaLeaf mini-app, which solves the heat equation in 2D using a matrix-free finite difference discretization [24].

Modifications to the classical EbE decomposition have recently been explored by Martínez-Frutos et al. [22]. The authors took a finer-grained approach to consider each degree of freedom (DoF) rather than a complete element, to formulate the DoF-by-DoF (DbD) method. It was shown in the context of elasticity problems that synchronization overhead from graph coloring is more costly than the unfavorable memory access patterns which it prevents, especially for 3D problems. Martínez-Frutos and Herrero-Pérez further explored the DbD method for elasticity problems on a fixed grid mesh so that only one local finite element matrix is required [20]. By varying material coefficients at the element level, between a constant “inside” and zero “outside”, different domain geometries were enforced on the same mesh. More recent work by the same authors applied these strategies to problems of robust topology optimization [21]. By leveraging the ability to simulate many different domain geometries, an optimal structural design for a given set of loading conditions can be found. Heterogeneous computing methods for topology optimization have also been studied for problems of designing domains with desired thermal properties. These studies have shown good results for problems involving steady state heat transfer in 2D [30] and 3D [21], using FEM solutions to the elliptic time-independent heat equation. Lastly, Martínez-Frutos and Herrero-Pérez have shown that multiple GPUs can be used effectively in the solution of topology optimization problems through task-level parallelism—the simultaneous evaluation of independent models that arise within a collocation strategy. All of the heterogeneous computing research described above, with the exception of TeaLeaf, was done with the CUDA platform.

The extension of heterogeneous computing beyond a single GPU has been an active research area, as further advances in scientific computing will require these paths. For example, Yang et al. study a general approach to partitioning the task of sparse matrix-vector multiplication across a range of heterogeneous architectures [31]. A strength of the approach is in emphasizing the mathematical model separate from the computational architecture, to allow greater flexibility. Similarly, Gao et al. optimize a matrix-explicit conjugate gradient solver for multiple GPU, with attention given to five candidate sparse matrix-vector multiplication formats to automatically select the one that is most appropriate for the hardware [7]. Considering that this is, by far, the most computationally expensive stage of the algorithm, the attention to the choice of memory storage format is warranted.

The present paper describes the theory and implementation of three approaches to a matrix-free preconditioned conjugate gradient algorithm for simulating transient heat conduction through a heterogeneous medium. Two are guided from previous studies, while our third combines benefits from both. The implementations differ in the interpretations of the DbD decomposition, as well as varying the use of a fixed grid or a general mesh, and in considering specialized hardware capabilities of GPUs. The most advanced implementation uses coalesced transactions with global memory for improvements in hardware efficiency, and is modified through a domain decomposition to run across dual GPUs. The domain decomposition is done in a way that minimizes communication between the two devices, so that the additional computational power can be effectively deployed. All of the implementation are made within the OpenCL computing framework, which is non-proprietary and free to use on any platform [28]. Scripting is done with the PyOpenCL package, providing readability and convenience with virtually no sacrifice in performance of the OpenCL API [16]. The performant code is available for public use, distribution, and modification [18].

3 Problem description

We present the mathematical and computational context for assembly-free finite element methods with a focus on the parabolic time-dependent heat equation PDE.

3.1 FE Formulation I (PDE)

We wish to solve the heat equation in three dimensions with spatially dependent material coefficients. In strong form, the boundary value problem is

[TABLE]

where the relevant thermal properties are the material density, $\rho$, specific heat, $C$, and thermal conductivity, $k$, all of which are assumed to be constant with respect to temperature. Discretizing in time with a $\theta$-scheme [19], the spatio-temporal temperature profile $T(\vec{x},t)$ is reduced to a finite set of temperatures at regularly spaced time increments, $\{T^{(i)}(\vec{x})\}_{i\in\mathcal{I}}$, $\mathcal{I}=\{0,1,\dots,t_{f}/\Delta t\}$. The problem is then converted to weak form by multiplying the strong form with a test function $\phi(\vec{x})$ and integrating by parts to give the operators

[TABLE]

We presume $T^{(i-1)}$ to be known, and find $T^{(i)}$ so that $a(T^{(i)},\phi)=L(\phi)$ for all $\phi$ in some family of functions. If this family is composed of a finite set of basis functions $\{\phi_{m}\}_{m\in\{1,\dots,N\}}$, the finite element method can be used to solve for $T^{(i)}$ as a linear combination of them: $T^{(i)}=\sum_{m}U^{(i)}_{m}\phi_{m}$, where $\vec{U}$ is a vector of coefficients. The discrete boundary integral of $f$ is computed to give the vector $\vec{F}$. The finite element method involves the resulting matrices of pairwise integrals of basis functions

[TABLE]

Computation of these matrices is the process of assembly. With $\mathbf{M}$ and $\mathbf{K}$ available, solving the weak form of the heat equation for all $\phi\in\{\phi_{m}\}_{m\in\{1,\dots,N\}}$ is equivalent to solving the matrix equation at each time step

[TABLE]

for $\vec{U}^{(i)}$. This differs from the process of solving for a steady state heat flow in two ways. First, the FEM system matrix has a more complicated structure, rather than only involving the stiffness matrix $\mathbf{K}$. Second, the system must be solved at every time step to produce a transient solution. The time discretization process requires its own considerations for numerical accuracy and stability [19]. In this work we set $\theta=0.5$, corresponding to a Crank-Nicolson method.

3.2 FE Formulation II (Assembly-Free Methods)

The matrices $\mathbf{M}$ and $\mathbf{K}$ are typically constructed by summing local contributions from each element in the assembly process. A local assembly matrix for element $e$, with $D$ degrees of freedom, contains the pairwise inner products of all basis functions with support in element $e$,

[TABLE]

Material property coefficients are taken to be constant over each element, so that the elemental assembly matrices depend only on the geometry of the domain. If all of the finite elements are the same size and shape, a single elemental assembly matrix can be reused and the mesh is said to have a fixed grid (FG). The assembly operator, \mathop{\mathchoice{\vbox{\Large\kern 1.23997pt\hbox{\mathsf{A}}\kern 1.23997pt}}{\vbox{\large\kern 1.03331pt\hbox{\mathsf{A}}\kern 1.03331pt}}{\vbox{\normalsize\kern 0.86108pt\hbox{\mathsf{A}}\kern 0.86108pt}}{\vbox{\small\kern 0.77498pt\hbox{\mathsf{A}}\kern 0.77498pt}}}\displaylimits, over the index set of elements $\mathcal{E}$ denotes the process of constructing a full system matrix from its local contributions. For example,

[TABLE]

The assembly operator can also be applied to contributions within a single vector over each element to give the full vector. It is only a notational convenience to describe the mapping from local degrees of freedom to sums over global degrees of freedom. As such, the following are valid notation for the general expression $\mathbf{M}\vec{x}=\vec{y}$

[TABLE]

meaning that assembly is computed in terms of the degrees of freedom (over index set $\mathcal{N}$) in the “outer loop” with each of their elemental contributions computed separately. The first method is an EbE approach, similar to the standard method of assembling $\mathbf{M}$. The second is a DbD approach, in which the necessary vector dot products are viewed with finer granularity [22]. The freedom of interpretation of the assembly operator gives rise to the different strategies for parallel matrix-vector multiplication that are described in Section 4.

To simplify notation, let $\mathbf{A}=\left[\mathbf{M}+\theta\Delta t\mathbf{K}\right]$ and $\mathbf{L}=\left[\mathbf{M}-(1-\theta)\Delta t\mathbf{K}\right]$. Set $\vec{b}=\mathbf{L}\vec{U}^{i-1}+\vec{F}$ and $\mathbf{A}_{e}=(\rho C)_{e}\mathbf{M}_{e}+\theta\Delta tk_{e}\mathbf{K}_{e}$ for $e\in\mathcal{E}$. Then the problem of finding $\vec{U}^{i}$ at each time step is reduced to solving

[TABLE]

We note once again that the explicit computation and storage of $\mathbf{A}$ and $\mathbf{L}$ is not necessary if $\{\mathbf{M}_{e}\}_{e\in\mathcal{E}}$ and $\{\mathbf{K}_{e}\}_{e\in\mathcal{E}}$ are available.

Assembly-free methods are especially useful if the spatially dependent material properties are not known in advance, or if many simulations are to be done over the same domain with varying coefficients. The generation of the mesh geometry and the computation of elemental assembly matrices can be done in advance and stored. Then all of the remaining computations required for a matrix-vector multiplication are parallelizable. Figure 1 illustrates a 3D FG mesh with tetrahedral elements and three sets of material properties parameterized by shading. Each of these spatially varying functions for the material properties have a simple functional form, and can be efficiently employed with elementary assembly matrices that are computed and stored beforehand, to treat changing analysis contexts under material property variation.

3.3 Preconditioned Conjugate Gradient

The matrix $\mathbf{A}=\left[\mathbf{M}+\theta\Delta t\mathbf{K}\right]$ is large, sparse, symmetric, and positive definite. The system above is thus solvable with the preconditioned conjugate gradient (PCG) algorithm [27]. At each time step, let $\vec{x}=\vec{U}^{(i)}$ (to free superscripts and the index $i$), be set at some initial guess. There also must be provided $\mathbf{A},\vec{b}$, a residual tolerance, and a preconditioner matrix $\mathbf{P}$ for which $\mathbf{P}^{-1}\vec{x}$ is easily computable and $\mathbf{P}^{-1}\mathbf{A}$ is relatively well conditioned, such as the Jacobi preconditioner. The algorithm involves the following linear algebra operations: matrix-vector multiplication (MVM), diagonal inverse matrix-vector multiplication (DIMVM), vector-vector multiplication (VVM), and adding scalar multiples of vectors (AXPY). The dominating computation is the MVM in each iteration. All of the other computations are easily parallelizable. The parallelization and implementation of the MVM on GPUs for system matrices of heat conduction FEM problems is the focus of this work. The subsequent PCG algorithm is shown in Algorithm 1 with each step annotated by its type of linear algebra operation.

A Jacobi preconditioner is used for all of the methods described in this work [27]. The construction of a Jacobi preconditioner is a straightforward process that lends itself to element-wise parallel computation. We note that the development of preconditioners that can be computed on a GPU is an active area of research [20, 8], but is not a focus of the present study. A comparison of this strategy with a serial implementation having a stronger preconditioner is made in Section 5.3.

3.4 OpenCL Heterogeneous Computing Framework

Here, we summarize the framework of OpenCL computing to introduce common language that will be used in Section 4.

3.4.1 Computation hierarchy

With OpenCL, GPUs are programmed with kernels; small bits of C code that are sent in parallel to the individual cores [28]. At any given time, each of the many cores is acting as a work item, which is the most granular operating unit in the hierarchy. Work items are structured in work groups, with as few as one work item per work group. The user is responsible for defining the sizes and dimensions of the hierarchical structure, so that at the time of execution, each work item is provided with unique identifying information and the generic kernel code. The identifying information is:

• •

global_id ranging from 0 to the total number of work items in each dimension,

• •

local_id ranging from 0 to the number of work items in a work group, in each dimension,

• •

group_id ranging from 0 to the total number of work groups in each dimension.

The total number of work items, total number of work groups, and size of each work group is also available. The kernel explicitly tells each work item how to contextualize itself within the larger problem, to determine what data must be loaded from memory or computed privately. At the end of the kernel, results of the independent granular computations are written to memory.

3.4.2 Memory hierarchy

The use of memory on a GPU dictates programming strategies and the success or failure of an algorithm. At fully efficient throughput, a single core can execute several floating point operations per clock cycle [4, 5]. However, accessing data from the “slow” global memory location can take 400-600 clock cycles. This alone warrants special attention to the OpenCL memory hierarchy.

Data that is loaded onto the GPU, or stored as the output of work items, must be stored in global memory, which has space on the order of gigabytes. During execution, all work items have their own small amount of private memory, which is on-chip and not visible to any other work item. This is fast to access, and used for variables that can take different values across every work item. In between, there is local memory, which is also on-chip, and shared among a single work group. There are usually tens of kilobytes reserved for local memory for each work item. Access to local memory is roughly 100 times faster than accessing global memory, provided that work items within the work group aren’t making conflicting calls (“bank conflict”). Local memory is allocated outside of the kernel, and cannot be freely initialized with specified data. Finally, kernels can consider certain data as constant memory. This data is physically still in global memory, but a kernel cannot write to it. When a kernel reads data from constant memory, it is cached, so that subsequent reads are fast. Global memory is not cached.

There are two strategies for gracefully managing reads from global memory when it is necessary. First, the latency can effectively be hidden if there is enough non-dependent computation to keep a work item busy between the time when the data are called and the time they are used. This is preferable, though not always possible. Second, a kernel can take advantage of the way that hardware loads data from global memory to private memory. A single transaction with global memory yields 32 words (such as 8 byte double-precision floats) of data, whether it is all called for or not. If work items that are indexed sequentially by global ID, request data from global memory that is organized in the same sequential way, the calls are automatically bundled and processed as one transaction. This process is the simplest form of coalesced memory access. There have been advances in hardware, and in OpenCL standards, to provide more flexibility, such as allowing permutations of the 32 sequential words to 32 work items, that are blocked together but not necessarily in the same order.

3.4.3 Programming strategy

The OpenCL programming approach is as follows:

Investigate the hardware to find and define a host (CPU/hard drive etc., the conventional computing environment), its devices(s) (GPU with its onboard memory), and define a context the overall computing environment. 2. 2.

Define initial variables on the host. 3. 3.

Load data onto a device. This includes reserving space for data that a kernel will write later and anything that is meant to last from one kernel to another. The benefit of defining all of the memory space in advance is that the user may specify whether each buffer is read or write only (or both) for both the host and device, or even if it is known that the host will never try to read it. Then the space that is allocated will be optimal for however the data will be treated. 4. 4.

Build the compute kernels. This involves compiling the C code and specifying pointers to memory buffers where its arguments can be found. A single kernel program can be built multiple times with different arguments, as is the case with the vector-vector multiplications in steps 5, 8, and 15 in Algorithm 1. Each of these are built independently. 5. 5.

Define a queue in the context. The host can enqueue kernels, memory transfer operations, or wait fences. OpenCL turns these into individual tasks that are performed as cores on the GPU become available. Code is written as if the context has infinitely many cores to run in parallel, and then the queue manages the execution of code on available hardware. Flags can be used with enqueued commands to ensure that all tasks from one kernel are finished before any tasks from the next kernel start, in case memory is being written and then read in a dependent way. 6. 6.

Enqueue commands to copy memory from a device to the host. This can be the final result of computations, or in the case of this work, the PCG residual, so that the host can decide whether or not to begin another iteration of enqueueing kernel commands.

4 Implementation of Assembly-Free Methods

We outline three assembly-free algorithms for matrix-vector multiplications. The differences arise from the flexibility in interpreting the assembly operator demonstrated in Equation (1). For each method, the explicit assembly equation is provided, along with an outline of the memory and computational hierarchy for a parallel implementation on GPUs. Further details are discussed in A so that broader concepts behind the implementations can be the present focus. We begin by giving context of the particular geometry of the problem.

4.1 Mesh Geometry

A 3D regular mesh with linear tetrahedral elements is generated based on the domain boundaries and the number of divisions in each dimension. The domain is then divided into rectangular prisms according to these divisions. The examples here are all cubes for simplicity. Each cube is then subdivided into six tetrahedra, as shown in Figure 1. A benefit of using tetrahedral elements is that their elemental assembly matrices have size $4\times 4$. The OpenCL specification allows dot products of the 4 element floating point vector data type, double4, with a single instruction. This is the central operation for all of our fast MVM methods, regardless of the interpretation of the assembly operator.

The six tetrahedra within each cube are indexed in a sequential way—six in the first cube, then six in the next cube in the $x$ direction, etc. until the $y$ dimension is incremented, and then the next slice in the $z$ dimension begins after that. An emphasized view of the six tetrahedra is given in Figure 2. Each cube has consistent indexing as to which of its eight corners correspond to each element within; allowing for a single reference table to be stored in memory. The structure of this kind of mesh also permits local computation, such as determining spatially dependent material properties, without the need for a table of $xyz$-locations of each vertex. This information can be determined from the one-dimensional index of the vertex along with a small table of the domain boundaries and the number of divisions in each dimension. Substituting small computation from general domain information, in place of loading the same data from a lookup table, is important for parallelizing an algorithm for a GPU.

4.2 Previous Strategies

We begin by providing the details of two approaches that use ideas from previous literature.

4.2.1 Implementation 1: Flexible DbD

The first implementation follows a flexible DbD strategy that does not assume that the FEM mesh lies on a fixed grid, but one which can be generated by a nonlinear transformation of one. Each global DoF of the vector $\vec{y}=\mathbf{A}\vec{x}$ in a general matrix equation is computed as

[TABLE]

where each pair of square brackets denotes a computation done by one work item in the implementation. Every elemental assembly matrix is stored in global memory, with a preprocessing step that determines $\{\mathbf{A}_{e}\}_{e\in\mathcal{E}}$ and $\{\mathbf{L}_{e}\}_{e\in\mathcal{E}}$ based on local material properties. A subset of this computation furnishes the Jacobi preconditioner matrix at the same time, by only storing diagonal entries. Work groups are responsible for each set of six elements in a cube, with 24 work items per work group, each corresponding to one element-DoF pair. Elemental assembly matrix coefficients are computed one time, and the scaled matrices are stored in double4 vectors in “element order” within global memory. In the computation, each of the 24 work items loads vertex data to local memory, including duplicated vertices, to alleviate conflicting simultaneous memory access from multiple work items. Then it reads its row from the elemental assembly matrix data, takes necessary entries from local vector data, performs a double4 dot product, and stores the result as its contribution to the global vertex in “vertex order”. In a second pass, one work group for each global vertex reads these 24 consecutive contribution and sums them into the single entry of the result vector. A similar two step approach is described by Martínez-Frutos et al. as an alternative to atomic incrementation to a single memory location for each degree of freedom [22].

The necessary data and computational responsibilities of the first pass of this process are illustrated in Figure 3. The left panel shows the scope of each work group—six elemental assembly matrices and eight entries of the input vector $\vec{x}$, which are loaded into the work group’s local memory. The right panel illustrates the responsibilities of each work item—one element-DoF pair, corresponding to a length-four vector dot product that it must compute and store. For clarity and consistency with the following sections, the necessary data for this single work item is emphasized in the left panel as well.

Splitting the assembly operator into two explicit sums requires a sort of data transpose at some point of the computation. This is because assembly matrix data is naturally stored in element order, while the output must be combined and stored in vertex order. The two step approach here ensures efficient data retrieved from global memory for both steps. The consequences of writing data to potentially distant locations at the end of the first step are hidden, since the data is enqueued to be written, and then the kernel can be restarted and the computation continues on while the writing takes place. Furthermore, we note the second step can be slightly modified to perform an extra vector-vector addition of the form $\vec{y}=\mathbf{A}\vec{x}+\vec{b}$ in the same kernel, which is one of the operations of the PCG algorithm.

4.2.2 Implementation 2: Single Pass FG DbD

The second implementation follows the coarser parsing of the assembly operator, with

[TABLE]

Once again, the square brackets denote the computations of a single work item, so that this method only requires a single pass. Single pass, fixed grid DbD strategies have been previously explored in References [20, 23]. The trade off here is that this requires more data to be loaded into local memory for each work group for the full determination of a degree of freedom. Consequently, more memory must be loaded overall. This is partially alleviated by setting the work groups to be as large as allowed by hardware limitations, so that most of the data are reused by adjacent degrees of freedom. Previous strategies to maximize on-chip memory have organized the input data into 2D square patches, which can be loaded from global memory in a coalescent way [26]. However, in order to have access to the data in neighboring elements, a halo of data around the patches must be loaded inefficiently. We choose to organize the necessary input data into 3 $\times$ 3 rectangular blocks, and as long as possible, so that all global memory access can be coalesced in the long direction, as shown in Figure 4. For our hardware, specified in Section 5, this corresponds to 64 work items and 27 KB of local memory usage. We also implement this approach to be used on a fixed grid mesh so that only one elemental assembly matrix is required. This allows the elemental assembly matrix to be stored in constant memory, so that it does not need to be loaded anew by each work item. The Jacobi preconditioner is computed in the same way as was described in Section 4.2.1, although the local material properties are not precomputed. This is not necessarily detrimental, since computing local material scaling coefficients is a small computation which can hide the latency of loading data from the input vector.

A visual overview of this method is provided in Figure 4. Each work item is responsible for all contributions to one entry in the output vector, so it needs data from 24 local assembly matrices and 27 entries of the input vector. An internal loop scans through all 24 element-DoF contributions, computing elemental scaling coefficients, performing dot products, and cumulatively summing the result.

4.3 Implementation 3: FG DbD with Memory Coalescing

The third implementation combines advantages from the first two. It is similar in structure to the first, except with increased responsibility to each work item, larger work groups, and the restriction to a fixed grid. The explicit interpretation of the assembly operator for this method is the same as Equation (2), although the full elemental assembly matrix-elemental vector multiplication is carried out by a single work item here. The primary characteristic of this third method is that work groups are structured so that all data that is loaded from global memory in the first pass is done so with a coalesced access pattern, as described in Section 3.4.2. In this sense, it is similar to the second implementation, except that the work group size is determined by the amount of data that can be loaded from a single coalesced memory read rather than by maximum local memory capacity. Furthermore, only four sections of memory are required for the first pass, as shown in Figure 5. The element-DoF contributions to each entry of the output vector are collected and summed as much as possible before writing back to global memory. While a second pass is still required to add the contributions among work groups, this process is faster than in the first implementation since there are fewer terms to sum. Additional details behind the efficient use of coalesced memory access now follow.

Since a fixed grid is assumed, the only data that must be loaded from global memory is the input vector. To do this efficiently, it must be loaded in blocks of 32 consecutive entries by blocks of 32 consecutive work items. A single memory read of this form does not give enough data to perform assembly for any element. However, if four blocks of memory are read, corresponding to the four horizontal edges of a long rectangular prism, then every tetrahedral element within that prism can be integrated over. That is, the contributions from each of these elements to the degrees of freedom that are loaded can be computed.

Using the memory access length of 32 as the guiding limit, 31 full cubes of six tetrahedral elements will be integrable. Therefore, work groups of 186 work items are invoked—one work item for each tetrahedral element. Since only 128 work items are necessary to read data from global memory, some work items are assigned a dual purpose and some remain idle during the loading process. Those first 128 receive both a global elemental index for their integration responsibility as well as a global vertex index for loading data from global memory into local memory. Figure 6 depicts the way in which global data is accessed by each work item for one of the four edges of the long rectangular prism. The other edges are treated in the same way, with offset information determined in-kernel based on how many divisions are made in the domain in each dimension. The contribution to a vertex can only be computed if all four vertices of an adjacent element are included within a work group. Therefore, only the 30 internal vertices from each block of 32 is contributed to in storage (except for the first work group, and possibly the last). Vertices on either end comprise a halo region that must be loaded for the computation within the local memory geometry. This partitioning pattern is also demonstrated in Figure 6. The special treatment of halo data on the region boundary has been explored through several GPU computing methods, both the context of the global domain and for local computation [24, 23, 26]. We next discuss our global memory halo.

Another consequence of coalesced memory access is that some consecutive elements of the input vector do not actually share an element in the domain. This is demonstrated in Figure 7 for a small example mesh. We avoid this problem by padding the domain with non-physical elements in the $+x$ and $+y$ directions from the perspectives of loading and computation, and then the contributions of these elements are ignored upon storage of the result. This produces some amount of wasted computation, but the relative volume of padding elements to the total mesh volume decreases as the granularity of the mesh increases. Furthermore, efficiency gains from coalesced memory access strongly outweigh the losses from this wasted computation.

4.4 FG DbD with Memory Coalescing on multiple GPUs

The fixed grid coalesced method described above is modified for use on dual GPU. The domain is split in the $z$ direction according to the additive Schwartz method [2]. This keeps all vertices in each subdomain in adjacent blocks of memory. A fraction of the domain to be assigned to device 1, $m$, is specified, such as 0.5. This fraction of $z$ slices, rounded up, with one additional layer, is the number of vertices that device 1 receives for computing, $m_{1}$. The rest of the vertices, in addition to two overlapping layers are assigned to device 2, totaling $m_{2}$. This split is demonstrated in Figure 8. If input vectors are initialized from a full set of global data, then one matrix-vector multiplication can be performed on each device, and the result can be faithfully reconstructed. For more than one sequential matrix-vector multiplication, the shared boundary data must be updated. Only one layer of vertices needs to be transferred in each direction. In the PCG algorithm, the solution vector $x$ is initialized at the beginning, and the intermediate vector $d$ must be transferred at each iteration. This memory transfer presents a bottleneck in the method, so that gains in speed are expected only with large systems for which MVM takes much longer than a GPU-to-GPU memory copy.

The scalar results of dot products must also be communicated between devices, both the residual $\delta$ and step size $\alpha$. To do this, each partial dot product is handled separately, omitting the extra boundary vertices. The partial results are stored in their own buffers, which are then transferred both directions. We thus require four memory buffers for each scalar quantity, one native buffer for the partial results on the device which computed it, and one target buffer on each device to store the copied value from the other. Many permutations of methods for the bidirectional communication of partial dot products were considered, such as copying the contributions to the host for summation. This method was found to be the fastest reliable way for each device to receive the full results. No special kernel is needed to combine them. Rather, modifications to the AXPY kernel are made to accept both buffers on a device and sum them as part of the existing process.

The only remaining matter for modification is adjusting the memory objects that contain vertex-to-spatial-location information. One each is made for the two devices, with the second encoding the $z$ shift for its first index. Then each device runs PCG in parallel, and can contextualize the location of vertices within the total domain. The host initialized two sets of every kernel that has been described for the serial method, as well as launching one queue to transfer boundary and scalar information and properly wait. Lastly, we note that the implementation here is for dual GPUs, but could be easily extended for increased distribution: with multiple GPUs. Additionally, the additive Schwartz method does not require uniformly sized decomposition regions. Therefore, heterogeneous GPUs could be used together, each taking a share of the total domain in order to balance the load across the devices.

5 Experiments and Discussion

The results shown here have been produced on AMD FirePro D700 GPU with 2048 streaming processors, 6GB of onboard memory, and up to 32KB of local memory per work group. The serial computations are done on a 2.7 GHz Intel Xeon E5 processor using a single core. The CPU algorithm is implemented in Python and optimized for solving heat conduction problems in many of the same ways as the GPU methods discussed previously. Elemental assembly matrices are precomputed so that assembly of the sparse system matrix can be done efficiently after the specification of material parameters, analogously to the parallel GPU algorithms. Furthermore, Jacobi preconditioning is specified so that the comparisons below are as fair as possible.

5.1 Performance Comparison

We first report the performance of each implementation over a range of problem sizes by simulating uniform heating on the front face of a two-layer laminate. The domain for this problem is a rectangular prism $\Omega=[-15,15]\times[-15,15]\times[0,10]$ with $f(\vec{x})=1$ on $x_{3}=0$ and zero everywhere else, and $T_{\text{ambient}}=0$. The material parameters are specified as

[TABLE]

corresponding to mild carbon steel and its solid corrosion products, assumed to be iron (III) oxide, Fe2O3. We perform 50 PCG solutions with $\Delta t=0.01$ seconds and a relative residual tolerance of $10^{-6}$. We vary the density of the tetrahedral mesh over the rectangular prism and measure the wall clock time per PCG iteration for each method. This provides a performance metric that accounts for computation time as well as the transfer of memory, and is consistent with the reporting of similar experiments [22, 20]. The results of these experiments are shown in Figure 9. Although there is sublinear progression with coarse meshes as memory transfer time dominates, all of the implementations exhibit linear increase in computation with the number of degrees of freedom. The rate of increase of wall clock time for each GPU implementation, taken from the linear asymptote, are compared with the sparse CPU method in Table 1. We note that these factors may differ from the factors taken from comparisons using similar implementations in the literature [22] for a number of reasons: lack of access to specific algorithmic details, hardware and API differences, and because our CPU baseline implementation is already highly optimized for this class of problem.

We also include, for comparison, the same profiling results when the GPU methods are set to compute with single-precision floating point arithmetic. As seen in Figure 9, the total timer per iteration decreases for all methods. We note that the point at which the dual GPU implementation becomes feasible is delayed to roughly double the problem size, as the balance between computation and memory transfer time is shifted. As the necessary computations become faster with single-precision arithmetic, the memory transfer overhead is more significant of a factor.

The experimentally observed linear scaling in computation with the number of degrees of freedom is expected based on the algorithms for computing MVM. Furthermore, it is seen that the efficiency of an implementation is directly related to the its care in the treatment of loading data from global memory. The second implementation has advantages over the first in that it does not require loading elementary matrix data or a second pass to sum contributions to a MVM. However, the cost is that it involves loading nine separate strips of data from the input vector to provide full information about all elements that are adjacent to a vertex. We see that this overhead outweighs the other merits of the implementation. The third implementation successfully adapts the memory management benefits of the previous two. A MVM requires only four strips of data from the input vector for the first pass, which are coalesced into four reads from global memory, and the contributions are partially summed so that the second pass will have less work. For the smallest problems considered here, the sparse matrix method is fastest, and the implementation on dual GPUs is the slowest. Further investigation into the dual GPUs implementation is done in the following section.

5.2 Multiple GPUs

Performing the third PCG method on dual GPUs gives better performance for large problems, as expected. We see in Figure 9 that for problems with fewer than about 400,000 degrees of freedom, the increased memory transfer costs dominate total computation time. For larger problem, the burden of computing MVMs grows, so that increased parallelism is worth these costs. We use OpenCL profiling operations to confirm that the time spent actively computing MVMs is reduced by the same factor as the workload sharing (including the overlapping region), even though the directly observed total time per iteration exhibits smaller reduction. The effective doubling of performance reported in Table 1 is taken from the relative slopes of the linear asymptote of each profiling curve. In other words, for the limit of problems in which computation dominates communication overhead. For problems larger than those considered here, it is expected that even more GPUs become practical, and we note that the extension is straightforward with the additive Schwartz method for domain decomposition.

5.3 Further CPU Comparison

We continue with a deeper comparison between our best performing heterogeneous computing method with a more sophisticated serial sparse matrix method. Since a standard CPU sparse matrix method can access the system matrix quickly, it is feasible to compute a better preconditioner matrix than the Jacobi preconditioner. This can take more time to produce and compute inverse matrix-vector multiplications with, but can vastly reduce the number of iterations required. We find that an incomplete Cholesky factorization with drop tolerance of $10^{-3}$ performs well for the simulations discussed in this paper [9]. Table 2 shows the benefits of using this preconditioner for a transient heat conduction boundary value problem over 50 time steps. Since it is no longer fair to compare the average time per PCG iteration, we report the total time that both methods take to provide a solution, starting from knowledge of the elemental assembly matrices and a parameterization of material properties, including the incomplete Cholesky factorization for the serial method. The effective time per iteration is determined from this adjusted total time, rather than exclusively from the time in the PCG outer loop. A graphical breakdown of the computations involved is presented in Figure 10.

We see that the sparse matrix CPU method with incomplete Cholesky preconditioning takes similar overall time for modestly large systems. However, the differences in wall clock time diverge as explicit storage of the system matrix becomes more demanding. The heterogeneous computing method developed in this research outperform even sophisticated serial algorithms.

5.4 3D Coefficient Inverse Problem

An important benefit of the currently proposed methods is that many simulations of heat conduction can be rapidly performed, in sequence, over a domain with varying thermal properties. This scenario arises in the solution of coefficient inverse problems, such as determining internal properties of a structure from noisy temperature measurements at its surface, in response to a known energy input. We explore such a case that has been earlier studied in two dimensions, motivated by a real-world corrosion detection problem [6]. By extending previous analysis methods from the literature to a full 3D model, we not only approach a more physically realistic use case, but are also able to relax certain symmetry restrictions on the boundary heat flux $f$. As a result, we found that an internal corrosion profile can be recovered with higher confidence, while using a heat source only one tenth as powerful as what was required in earlier work [6]. The details of the numerical experiment are now discussed.

Corrosion may form in a bridge structure in the crevice where a lower truss chord member meets a flat steel gusset plate connection element. It is observed that this corrosion within the gusset plate has a well defined geometric form, constant along the horizontal length of the truss chord, and varying as a quadratic function in the vertical direction [6]. The severity of the corrosion is parameterized by the penetration depth of the apex of the parabola into the steel, $\theta$, illustrated in Figure 11. We consider a case of corrosion that has penetrated 3.175 mm into a plate that is 12.7 mm thick: corresponding to a 25% section loss. The relevant thermal properties are the same as in Section 5.1: A thermal input is furnished in the form of a 10 W laser beam having a Gaussian profile with 2 mm beam width. The structure is heated for $T_{F}=10$ seconds, after which the surface temperature is recorded with an thermal camera. Surrogate field data, $\mathcal{D}$, are created with a high fidelity forward model (1.16 million degrees of freedom and 0.01 second time steps), interpolating the FEM solution to a finer rectangular grid, averaging the temperatures over each pixel area, and contaminating the result with noise; all so as to approximate a plausible digital thermal camera measuring device. For consistency with previous work, the camera is assumed to follow a noise model consisting of contaminating each pixel with independent and identically distributed Gaussian noise with zero mean and standard deviation of 0.1 ∘C, then rounding the result to the nearest 0.1 ∘C.

The inverse problem is solved with Markov chain Monte Carlo (MCMC), a Bayesian inference method [10]. The solution comes in the form of dependent samples taken from a probability distribution over the corrosion penetration depth, $p(\theta|\mathcal{D})$, known as the posterior distribution. Since there is randomness in the surface temperature measurements, MCMC is able to automatically propagate a measure of uncertainty to the posterior distribution. The mean and standard deviation of a large number of posterior samples provide useful information about the underlying corrosion depth. The consequence is that every sample comes at the cost of performing a simulation with the candidate corrosion depth, $\hat{\theta}$, and evaluating the probability that the resulting FEM solution gave rise to the observed temperature data after passing through the camera noise model. This is known as the likelihood, $p(\mathcal{D}|\hat{\theta})$. Assuming an “uninformative” uniform distribution over $\theta$, that is called a prior, the corrosion depth is equally likely to be anywhere within the gusset plate thickness, Bayes’ theorem states that the posterior distribution is proportional to the likelihood distribution, $p(\theta|\mathcal{D})\propto p(\mathcal{D}|\theta)$, only scaled by a constant. Thus, the problem of estimating the posterior distribution can be reduced to sampling candidate values for $\theta$ and evaluating their corresponding likelihoods.

Markov chain Monte Carlo provides a systematic algorithm for producing the desired samples over $\theta$. Starting from an initial guess, $\theta_{0}$, a new candidate, $\hat{\theta}$, is randomly generated from some small neighborhood of $\theta_{0}$. The likelihoods for each value, given the observed data, are calculated through the use of the FEM solver. If the new sample is more likely to have produced the data, then it is set as $\theta_{1}$. Even if $\hat{\theta}$ was less likely to have produced the observed data, it may be randomly selected as the next sample with probability $p(\mathcal{D}|\hat{\theta})/p(\mathcal{D}|\theta_{0})$. Otherwise, $\theta_{0}$ advances, repeated as sample $\theta_{1}$. This algorithm, known as Metropolis-Hastings sampling [10], is summarized as

[TABLE]

Under some technical conditions that are easily verified for our system [17], the sequence $\{\theta_{i}\}_{i=1,2,3,...}$ is guaranteed to converge to samples from the desired distribution $p(\theta|\mathcal{D})$. The sampling process is summarized in Figure 12. The FEM implementation developed in this work is nested within two loops, as the critical stage of computation. Our GPU approach is massively parallelized for this task, and is designed to provide rapid successive solutions with only the transfer of $\theta_{i}$ from the host to the device (i.e. very low communication overhead). Furthermore, if the elemental assembly data is precomputed and loaded onto the device, the entire process can be carried out with minimal memory transfer.

It is customary when using MCMC to discard some number of initial samples, to allow the chain to “burn in” to the posterior distribution, and away from its arbitrary initial value. For these experiments, we perform 200 burn-in iterations, starting from the middle of the prior distribution, and then record the subsequent 2500 samples. The trajectory of the Markov chain and a histogram of the samples are shown in Figure 13. It can be seen that the samples fall tightly around the ground truth value of 3.175 mm, denoting a confident solution to the coefficient inverse problem. The mean estimate is $3.16$ mm with standard deviation $0.05$ mm.

The demonstrated sequence of transient heat conduction simulations completed in 9 hours. If the same simulations are attempted using FEniCS, a modern open-source FEM computing platform [19], the system runs out of memory in the assembly process. Computing on a series of smaller problems and extrapolating the run times (ignoring all contributions from just-in-time compilation) yields an expected total time of 142 hours. Using the in-house sparse matrix CPU code that was written specifically to solve this kind of sequence of heat conduction simulation, and discussed in Section 5.3, would have taken approximately 38 hours.

6 Conclusions

In this paper we described the development over three implementations of fine-grained assembly-free finite element simulation of heat transfer through heterogeneous media with the use of heterogeneous computing. The methods are motivated in terms of the interpretation of the FEM assembly operator as well as their particular implementation on GPUs. The implementation which gave the best performance was adapted for its simultaneous application on multiple GPUs. The ease of application of these methods to problems with spatially varying material properties is a direct consequence of the decomposition of the problem into local geometric components, rather than focusing on a global system matrix. Furthermore, successive simulations have low overhead, as mesh generation and computation of elemental assembly matrices can be done in advance. For these reasons, the methods that are developed in this research are particularly well-suited as forward models for coefficient inverse problems and topology optimization problems involving transient heat conduction.

We finally note that the algorithm design methods that we have described are not restricted to the heat equation, and can be used as guidelines for solving other PDEs with spatially varying coefficients. To promote the use and extension of our work, we have used the non-proprietary OpenCL API and made the code for our best-performing implementation (both single and dual GPU versions) publicly available [18].

7 Acknowledgments

This work was supported by the National Science Foundation (DGE-1144153).

Appendix A GPU PCG Algorithm Details

This appendix contains descriptions of the custom kernels that were written to perform linear algebra operations, how they use GPU memory, and how they fit into the PCG algorithm.

A.1 Memory

We first introduce the FEM data and indexing variables that are used throughout our kernels:

• •

M and F: arrays of data for elemental assembly matrices. Either the single matrices for fixed grid methods, or all entries for every element for the general methods

• •

DoFMapLocal and coordinateMapLocal: small arrays of indices based on the geometry of the mesh to quickly determine local DoF to element assignment

• •

C: array of the number of divisions of the domain in each dimension

• •

vert_scale: array of the minimum vertex position and spacing between vertices in each dimension. Along with C, this allows the determination of the absolute location of a vertex given its global index. For a non-fixed grid approach, the user may start with a uniform mesh and deform it in a predetermined way (e.g. quadratically scaling the vertex locations), so that information can be used inside a kernel to still enable the recovery of absolute vertex position.

• •

corr_bound: spatial information about what part of the domain is corroded (or in general, a different material). For uniform corrosion after a certain depth, this can be that distance in the $z$-direction. For elliptical corrosion pits, it could contain the coordinates of the center of the ellipse and its axis lengths. The kernel must be programmed to know how to interpret this.

• •

mat_coefs: values for $\rho C$ and $k$ for the different materials

Memory is allocated in GPU global memory buffers with flags to specify how the host and the device will access them. The flags are self-explanatory, and are combined with logical OR. The combinations used here are:

• •

HOST_TO_DEVICE_COPY = (READ_ONLY $|$ HOST_WRITE_ONLY $|$ COPY_HOST_PTR)

• •

HOST_TO_DEVICE_USE = (READ_ONLY $|$ HOST_WRITE_ONLY $|$ USE_HOST_PTR)

• •

HOST_READ_WRITE = (READ_WRITE $|$ COPY_HOST_PTR)

• •

PINNED = (READ_WRITE $|$ USE_HOST_PTR)

• •

DEVICE_READ_WRITE = (READ_WRITE $|$ HOST_NO_ACCESS)

The memory buffers allocated for the FG DbD method are enumerated in Table 3. Other methods do not differ much at this level.

A.2 Kernels

A.2.1 VVM_A

First stage of a vector-vector multiplication. Takes pairwise scalar products of vector elements and then uses a standard parallel “reduction” algorithm to sum the results in $\mathcal{O}(\log(N))$ complexity. The work group size is maximized for the available hardware, and the partial sum is reduced by a factor of two at each step. The total number of elements in the partial sum can only be reduced by a factor of the maximum work group size with a single kernel call. In case the length of the vector, nVertices, is not a multiple of the maximum work group size, max_wg_size, the number of work groups is rounded up. The size of the result (number of elements) is then

[TABLE]

A.2.2 VVM_reduce

Intermediate and final stages of vector-vector multiplication. Takes an intermediate partial sum from VVM_A or itself and reduces it by a factor of max_workgroup_size. If nVertices $\leq$ max_workgroup_size, the vector-vector multiplication is completed. Otherwise further calls with this kernel are made. Since this can be iterative, set

[TABLE]

A.2.3 VVM_C

An alternative final stage for vector-vector multiplication for step 5 of PCG. Rather than storing the final result of the sum, the scalar delta_new is loaded and alpha = delta_new/(dTq) is stored, along with negative alpha.

A.2.4 VAVSP

Computes the elementwise sum of a vector and a scalar multiplied by another vector. The scalars are loaded from GPU global memory, so a negative scalar must be used if vector subtraction is desired.

A.2.5 DIMVM

Computes the matrix-vector multiplication where the matrix is the inverse of a diagonal matrix P. Element i of the result is x_i/P_{i,i}.

A.2.6 Beta_update

Computes the coefficient beta and updates delta in storage. Used to avoid data transfer between device and host for the performance of a small calculation.

A.2.7 u0_update

Computes an initial guess for the next time step, u_0+, based on a linear extrapolation from the initial guess of the current time step, u_0 and the PCG solution of the current time step, u_new,

[TABLE]

A.2.8 FGDbDMVM_A

The implementation of matrix-vector multiplication (including cases for $\mathbf{A}\vec{x}$ and $a\mathbf{A}\vec{x}+\vec{b}$) and the determination of $P$ is dependent on the assembly perspective. The kernels for the third implementation (FG DbD with memory coalescing) are discussed here.

The kernel computes contributions to a matrix-vector multiplication from each element according to coalesced memory access limitations. Data is loaded from the input vector as described in Section 4.3 and stored in local memory. Then each work item determines its global element index by first finding the cube in which it belongs: global_id/6 - group_id (global_id is an unsigned integer, so integer division automatically rounds down), and then the tetrahedron within the cube: local_id mod 6.

The cube index also corresponds to the global vertex index for its first corner. This is used with the arrays C and vert_scale to determine the vertex’s absolute position. Based on this, the array coordinateMapLocal helps give the $x,y$, and $z$ positions of the other vertices of the tetrahedron, so that they can be compared with corr_bounds to determine for each vertex whether it is in the corroded region of the domain or not. Material coefficients are taken from mat_coefs and averaged over the element for both $\mathbf{M}_{e}$ and $\mathbf{K}_{e}$.

Next, each element reads the data it needs about the input vector from local memory, referring to DoFMapLocal for the proper indices. This has been delayed as long as possible to hide the latency from the global memory load. Dot products are taken with the local assembly matrices, and the results are summed with the proper coefficients in another local memory array, Ax_split_local, which has 12 entries for every DoF, corresponding to the 12 possible contributing elements in the $\pm x$ directions. Finally, the work items that are responsible for loading and storing data sum over Ax_split_local for their DoF and write to global memory. The resulting memory buffer has 4 entries for every global DoF, one for every quadrant in the $y$-$z$ plane. These are all filled out in turn by further work groups.

The total number of work items needed for this kernel is found by determining the total number of elements that need to be considered including the padding in $+x$ and $+y$ directions, dividing by 30 since every work group yields the elemental contributions for blocks of 30 vertices, and multiplying by the work group size. As with vector-vector multiplications, we round up the integer division

[TABLE]

A.2.9 FGDbDMVM_B

Finishes the matrix-vector multiplication started by FGDbDMVM_A. Work items sum the four contributions to each DoF from Ax_split and store them in the final result array.

A.2.10 FGDbDMVM_C

Finishes the matrix-vector multiplication started by FGDbDMVM_A with an extra SAXPY operation so that a separate call to AXPY is not necessary. Work items sum the four contributions to each DoF from Ax_split, multiply them by a scalar, add them to an element from another vector, and store the result.

A.2.11 Jacobi_A

Determines contributions to the Jacobi preconditioner $P$. The algorithm is the same as FGDbDMVM_A, except instead of taking dot products with elemental assembly matrices and an input vector, the diagonal elements of the elemental assembly matices are scaled according to material coefficients, combined, and stored. Calling FGDbDMVM_B on the result produces the diagonal of $P$.

Arguments for this kernel are the same as for FGDbDMVM_A, except without the need for *x and *x_local.

A.3 PCG Again

The kernels and memory usage is as follows, folowing the steps in Section 3.3. The syntax below is a small abbreviation of the actual PyOpenCL code, and has the form

kernel_instance = kernel_name(args)

First compute $P$ and the right hand side vector $\vec{b}=\mathbf{L}\vec{u}_{i}+\vec{F}$.

• •

knl_PA = Jacobi_A(M_FG_buf, K_FG_buf, DoFMapLocal_buf, coordinateMapLocal_buf, C_buf, vert_scale_buf, corr_bounds_buf, mat_coefs_buf, theta, dt, Ax_split_local_buf, P_split_buf)

• •

knl_PB = FGDbDMVM_B(P_split_buf, P_buf)

• •

knl_RHS_A = FGDbDMVM_A(M_FG_buf, K_FG_buf, u0_buf, DoFMapLocal_buf, coordinateMapLocal_buf, C_buf, vert_scale_buf, corr_bounds_buf, mat_coefs_buf, (1-theta), dt, x_local_buf, Ax_split_local_buf, Ax_split_buf)

• •

knl_RHS_B = FGDbDMVM_C(Ax_split_buf, Fdt_buf, np.float64(1), b_buf)

Setup for PCG

• •

knl_1A = FGDbDMVM_A(M_FG_buf, K_FG_buf, x_buf, DoFMapLocal_buf, coordinateMapLocal_buf, C_buf, vert_scale_buf, corr_bounds_buf, mat_coefs_buf, theta, dt, x_local_buf, Ax_split_local_buf, Ax_split_buf)

• •

knl_1B = FGDbDMVM_C(Ax_split_buf, b_buf, np.float64(-1), r_buf)

• •

knl_2 = DIMVM(P_buf, r_buf, d_buf)

• •

knl_3A = VVM_A(r_buf, d_buf, nVertices, VVM_loc_buf, r1_buf)

• •

knl_3B = VVM_reduce(r1_buf, r1_size, VVM_loc_buf, r2_buf)

• •

knl_3C = VVM_reduce(r2_buf, r2_size, VVM_loc_buf, delta_buf)

Perform one PCG iteration

• •

knl_4A = FGDbDMVM_A(M_FG_buf, K_FG_buf, d_buf, DoFMapLocal_buf, coordinateMapLocal_buf, C_buf, vert_scale_buf, corr_bounds_buf, mat_coefs_buf, theta, dt, x_local_buf, Ax_split_local_buf, Ax_split_buf)

• •

knl_4B = FGDbDMVM_B(Ax_split_buf, q_buf)

• •

knl_5A = VVM_A(d_buf, q_buf, nVertices, VVM_loc_buf, r1_buf)

• •

knl_5B = VVM_reduce(r1_buf, r1_size, VVM_loc_buf, r2_buf)

• •

knl_5C = VVM_C(r2_buf, delta_buf, r2_size, VVM_loc_buf, alpha_buf, neg_alpha_buf)

• •

knl_6 = AXPY(x_buf, d_buf, alpha_buf, x_buf)

• •

knl_7A = FGDbDMVM_A(M_FG_buf, K_FG_buf, x_buf, DoFMapLocal_buf, coordinateMapLocal_buf, C_buf, vert_scale_buf, corr_bounds_buf, mat_coefs_buf, theta, dt, x_local_buf, Ax_split_local_buf, Ax_split_buf)

• •

knl_7B = FGDbDMVM_C(Ax_split_buf, b_buf, np.float64(-1), r_buf)

• •

knl_7 = AXPY(r_buf, q_buf, neg_alpha_buf, r_buf)

• •

knl_8 = DIMVM(P_buf, r_buf, s_buf)

• •

knl_9A = VVM_A(r_buf, s_buf, nVertices, VVM_loc_buf, r1_buf)

• •

knl_9B = VVM_reduce(r1_buf, r1_size, VVM_loc_buf, r2_buf)

• •

knl_9C = VVM_reduce(r2_buf, r2_size, VVM_loc_buf, delta_new_buf)

• •

knl_10 = Beta_update(delta_new_buf, delta_buf, beta_buf)

• •

knl_11 = AXPY(s_buf, d_buf, beta_buf, d_buf)

• •

knl_u0 = u0_update(u0_buf, x_buf)

The CPU handles logical decisions and calls these kernels according to the algorithm until convergence is realized.