# TAC+: Drastically Optimizing Error-Bounded Lossy Compression for 3D AMR
Simulations

Daoce Wang, Jesus Pulido, Pascal Grosset, Sian Jin, Jiannan Tian,
Kai Zhao, James Arhens, and Dingwen Tao  Daoce Wang, Sian Jin, Jiannan Tian,
and Dingwen Tao (corresponding author) are with Indiana University,
Bloomington, IN 47405, USA. Jesus Pulido, Pascal Grosset, and James Arhes are
with Los Alamos National Laboratory, Los Alamos, NM 87545, USA. Kai Zhao is
with University of Alabama at Birmingham, AL 35294, USA.

###### Abstract

Today’s scientific simulations require a significant reduction of data volume
because of extremely large amounts of data they produce and the limited I/O
bandwidth and storage space. Error-bounded lossy compression has been
considered one of the most effective solutions to the above problem. However,
little work has been done to improve error-bounded lossy compression for
Adaptive Mesh Refinement (AMR) simulation data. Unlike the previous work that
only leverages 1D compression, in this work, we propose an approach (TAC) to
leverage high-dimensional SZ compression for each refinement level of AMR
data. To remove the data redundancy across different levels, we propose
several pre-process strategies and adaptively use them based on the data
characteristics. We further optimize TAC to TAC+ by improving the lossless
encoding stage of SZ compression to efficiently handle many small AMR data
blocks after the pre-processing. Experiments on 8 AMR datasets from a real-
world large-scale AMR simulation demonstrate that TAC+ can improve the
compression ratio by up to 4.9$\times$ under the same data distortion,
compared to the state-of-the-art method. In addition, we leverage the
flexibility of our approach to tune the error bound for each level, which
achieves much lower data distortion on two application-specific metrics.

##  1 Introduction

The increase in supercomputer performance over the past decades has been
insufficient to solve many challenging modeling and simulation problems. For
example, the complexity of solving evolutionary partial differential equations
scales as $\Omega(n^{4})$, where $n$ is the number of mesh points per
dimension. Thus, the performance improvement of about three orders of
magnitudes over the past 30 years has meant just a 5.6$\times$ gain in spatio-
temporal resolution [1]. To address this issue, many high-performance
computing (HPC) simulation packages [2] (such as AMReX [3] and Athena++ [4])
use Adaptive Mesh Refinement (AMR)—which applies computation to selective
regions of most interest—to increase resolution. Compared to the method where
a high resolution is applied everywhere, the AMR method greatly reduces the
computational complexity and storage overhead; thus, it is one of the most
widely used frameworks for many HPC applications [5, 6, 7, 8].

Although AMR can save storage space to some extent, AMR applications running
on supercomputers still generate large amounts of data, bringing challenges to
data transmission and storage. For example, one Nyx simulation [9] with a
resolution of $4096^{3}$ (i.e., $0.5\times 2048^{3}$ mesh points in the coarse
level and $0.5\times 4096^{3}$ in the fine level ) can generate up to 1.8 TB
of data for a single snapshot; a total of 1.8 PB of disk storage is needed
assuming running the simulation 5 times with 200 snapshots dumped per
simulation. Therefore, reducing data size is necessary to lower the storage
overhead and I/O cost and improve the overall application performance for
running large-scale AMR simulations on supercomputers.

A straightforward way to address this issue is to use data compression.
However, traditional lossless compression techniques such as GZIP [10] and
Zstandard [11] can only provide a compression ratio by up to 2$\times$ for
scientific data [12]. On the other hand, a new generation of lossy compressors
that can provide strict error control (called “error-bounded” lossy
compression) has been developed, such as SZ [13, 14, 15], ZFP [16], MGARD
[17], and TTHRESH [18]. Using those error-bounded lossy compressors,
scientists can achieve relatively high compression ratios while minimizing the
quality loss of reconstructed data and post-analysis, as demonstrated in [19,
20, 21, 22, 23, 24, 25, 26].

Only a few existing contributions have investigated error-bounded lossy
compression for AMR applications and datasets. A common approach is to
generate uniform-resolution data by up-sampling the coarse-level data and
merging them with the finest-level data and then performing compression on the
merged data. However, this approach introduces redundant information to the
data, which significantly degrades the compression ratio, especially when the
up-sampling rate is high or there are multiple coarse levels to up-sample.

Recently, Luo et al. introduced zMesh [27], a technique that groups data
points that are mapped to the same or adjacent geometric coordinates such that
the dataset is smoother and more compressible. However, since zMesh maps data
points from different AMR levels to adjacent geometric coordinates and
generates a 1D array, it cannot adopt 3D compression which most HPC
simulations use. Moreover, zMesh is designed for patch-based AMR data111The
patch-based AMR data redundantly saves the data block to be refined at the
next finer level in the current coarse level (will be introduced in detail in
Section 2.3). with redundancy across different AMR levels to improve the
compression ratio. However, these coarser levels of redundant data are often
not used for post-analysis or visualization and hence can be directly
discarded to improve the compression ratio. For this case, the reorganization
approach proposed by zMesh cannot improve the data smoothness appropriately
(will be illustrated in Section 4).

To solve these issues, we propose TAC that removes the redundant data in
coarser level(s) and employs 3D lossy compression for each level. We note that
each level may contain many empty/zero regions, where data points are saved in
other levels, which may significantly decrease the data smoothness and hence
reduce the compression ratio. To this end, TAC either removes these empty
regions using adaptive partition strategy or partially pads them with
appropriate values, based on the density of empty regions. TAC also has an
optimization to reduce the time complexity of removing empty regions.

Although our partition strategies can remove empty regions, there are still
several challenges. Specifically, the partition strategies can generate many
(e.g., 3,000+) small data blocks with totally different shapes (e.g., 10
different shapes), whereas the SZ compressor performs poorly on relatively
small data sets. This is because the SZ compressor uses thousands of Huffman
trees to encode these small blocks separately, leading to a low encoding
efficiency. A naive solution to the heavy Huffman encoding cost is to
linearize/merge the thousands of small blocks into fewer larger blocks and
then pass these larger blocks to SZ. This approach can reduce the overhead of
the Huffman trees for encoding and hence increase the amount of data encoded
together.

However, TAC still has two limitations: (1) Most of the merged small blocks
are not adjacent in the original dataset, leading to rapid changes in the data
values at the boundaries of these non-neighboring blocks, which impacts the
accuracy of SZ’s predictor negatively. (2) Although we can compress data
blocks with the same shape together, the SZ compressor must be called multiple
times for each shape, resulting in the inevitable low performance of Huffman
encoding. To address these limitations, we further optimize TAC to TAC+ by
designing a Shared Huffman Encoding (SHE) approach for the SZ compressor. This
approach allows individual predictions for each small block while being
encoded using a single shared Huffman tree, which can improve the prediction
accuracy and compression ratio accordingly.

The main contributions are summarized as follows.

  * •

We propose to leverage 3D SZ compression to compress each level of an AMR
dataset separately. We propose a hybrid compression approach based on the
following three pre-process strategies and data characteristics.

  * •

We propose an optimized sparse tensor representation to efficiently partition
data and remove empty regions for sparse AMR data.

  * •

We propose an enhanced $k$-d tree approach to reduce the time overhead of
removing empty regions.

  * •

We propose a padding approach to improve the smoothness and compressibility of
dense AMR data.

  * •

We employ the SHE approach in the SZ compressor to reduce the high time and
storage costs of compressing multiple small blocks after the partition.

  * •

We tune the error bound for each AMR level to further improve the compression
quality in terms of two application-specific post-analysis metrics.

  * •

Experiments show that, compared to the state-of-the-art approach zMesh, our
proposed AMR compression can improve the compression ratio by up to
4.9$\times$ under the same data distortion on the tested datasets.

We evaluate our proposed compression method on eight datasets from three real-
world AMR simulation runs. The AMR simulations are well-known, open-source
cosmology simulations—Nyx [9]. We compare our method with four baselines
including zMesh using generic metrics such as compression ratio and peak
signal-to-noise ratio (PSNR) and application-specific metrics such as power
spectrum and halo finder. Our code is available at
https://github.com/FabioGrosso/3dAMRcomp.

In Section 2, we present background information about error-bounded lossy
compression, AMR method, $k$-d tree, and related work on AMR data compression.
In Section 3, we describe our proposed pre-process strategies, SHE approach,
and hybrid compression. In Section 4, we show the experimental results on
different AMR datasets. In Section 5, we conclude our work and discuss the
future work.

##  2 Background and Related Work

###  2.1 Lossy Compression for Scientific Data

There are two main categories for data compression: lossless and lossy
compression. Compared to lossless compression, lossy compression can offer a
much higher compression ratio by trading a little bit of accuracy. There are
some well-developed lossy compressors for images and videos such as JPEG [28]
and MPEG [29], but they do not have a good performance on the scientific data
because they are mainly designed for integers rather than floating points.

In recent years there is a new generation of lossy compressors that are
designed for scientific data, such as SZ [13, 14, 15], ZFP [16], MGARD [17],
and TTHRESH [18]. These lossy compressors provide parameters that allow users
to finely control the information loss introduced by lossy compression. Unlike
traditional lossy compressors such as JPEG [28] for images (in integers), SZ,
ZFP, MGARD, and TTHRESH are designed to compress floating-point data and can
provide a strict error-controlling scheme based on the user’s requirements.
Generally, lossy compressors provide multiple compression modes, such as
error-bounding mode and fixed-rate mode. Error-bounding mode requires users to
set an error type, such as the point-wise absolute error bound and point-wise
relative error bound, and an error bound level (e.g., $10^{-3}$). The
compressor ensures that the differences between the original data and the
reconstructed data do not exceed the user-set error bound level.

In this work, we focus on the SZ lossy compression (2021 R&D 100 Award Winner
[30]) because SZ typically provides a higher compression ratio than ZFP [22,
31] and higher (de)compression speeds than MGARD [31, 32] and TTHRESH [18]. SZ
is a prediction-based error-bounded lossy compressor for scientific data. It
has three main steps: (1) predict each data point’s value based on its
neighboring points by using an adaptive, best-fit prediction method; (2)
quantize the difference between the real value and predicted value based on
the user-set error bound; and (3) apply a customized Huffman coding and
lossless compression.

###  2.2 AMR Method and AMR data

Figure 1: Visualization (one zoom-in 2D slice) of three key timesteps
generated from an AMR-based cosmology simulation. The grid structure changes
with the universe’s evolution. The red boxes indicate different resolutions
within one AMR level.

AMR is a method of adapting the accuracy of a solution by using a non-uniform
grid to increase computational and storage savings while still achieving the
desired accuracy. AMR applications change the mesh or spatial resolution based
on the level of refinement needed by the simulation and use finer mesh in the
regions with more importance/interest and coarser mesh in the regions with
less importance/interest. Figure 1 shows that during an AMR run, the mesh will
be refined when the value meets the refinement criteria, e.g., refining a
block when its norm of the gradients or maximum value is larger than a
threshold.

Figure 2: A typical example of AMR data storage and usage.

Clearly, the data generated by an AMR application are hierarchical data with
different resolutions. The data of each AMR level are usually stored
separately (e.g., in a 1D array). For example, Figure 2 (left) shows a simple
example of two-level AMR data; “0” means high resolution (the fine level) and
“1” for low resolution (the coarse level). When the AMR data are needed for
post analysis or visualization, users will typically convert the data from
different levels to a uniform resolution. In the previous example, we will up-
sample the data in the coarse level and combine it with the data in the fine
level, as shown in Figure 2 (right).

###  2.3 Tree-based and Patch-based AMR Data

There are two types of techniques to represent AMR data: patch-based AMR and
tree-based AMR [33]. The main difference between them is that the patch-based
AMR technique generates AMR data with redundancy across different levels. In
other words, the patch-based AMR data structure redundantly saves data blocks
to be refined at the next level in the current level, simplifying the
computation in the refinement process. By comparison, the tree-based AMR
technique organizes the grids on the tree leaves, so there is no redundant
data across different levels. But tree-based AMR data is more complex for post
analysis and visualization compared to patch-based AMR data [34].

In this work, we focus on a state-of-the-art patch-based AMR framework AMReX.
Note that since the redundant coarser-level data in the patch-based AMR will
not often be used in post-analysis, we discard them during compression to
improve the compression ratio.

###  2.4 Existing AMR Data Compression

1D AMR Compression: The main challenge for AMR data compression is that the
AMR data is comprehensive and hierarchical with different resolutions. A naive
approach is to compress the 1D data of each AMR level separately. However,
this approach loses most of the topological/spatial information, which is
critical for data compression. zMesh [27] is a state-of-the-art AMR data
compression based on the 1D approach. Different from the naive 1D approach,
zMesh re-organizes the 1D data based on each point’s coordinate in the 2D
layout; in other words, zMesh puts the points neighbored in the 2D layout
closer in the 1D array. It can increase the data smoothness/compressibility to
benefit the following 1D compression such as SZ on patch-based AMR data with
redundancy across different AMR levels. However, zMesh does not leverage high-
dimensional compression, while many previous studies [14, 35] proved that
leveraging more dimensional information (e.g., spatial/temporal information)
can significantly improve the compression performance. Moreover, it only
focuses on 2D AMR data. Our work aims to leverage high-dimensional data
compression and supports 3D AMR data.

High-dimensional AMR Compression: Similar to the idea described in Section
2.2, a straightforward way to leverage 3D compression on 3D AMR data is to
compress different levels together by up-sampling coarse levels. However, this
approach must handle extra redundant data generated by the up-sampling
process. As shown in Figure 2, 1A, 1B, and 1C are redundant points in the
compression. Note that the storage overhead of these redundant points will be
higher when more data are in the coarse levels or the up-sampling rate is
higher, especially for 3D AMR data. This is because we only need to duplicate
one point from the coarse level 4 times for 2D AMR data but 8 times for 3D AMR
data, with an up-sampling rate of 2. Another limitation of this approach is
that it cannot apply different compression configurations (e.g., error bound)
to different AMR levels. This is because after up-sampling all data points
will have the same importance. However, the purpose of using the AMR method is
to set different interests to different AMR levels, so the error bound for
each AMR level can be chosen adaptively.

###  2.5 $k$-D Tree for Particle Data Compression

$k$-d tree [36] is a binary tree in which every node represents a certain
space. Without loss of generality, for the 3D case, every non-leaf node in a
$k$-d tree splits the space into two parts by a 2D plane associated with one
of the three dimensions. The left subspace is associated with the left child
of the node, while the right subspace is associated with the right child.
$k$-d tree is commonly used in particle data compression [37, 38, 39] to
locate each particle and remove empty regions. Specifically, a $k$-d tree
keeps dividing the space in between along one dimension until the space is
empty or contains only one particle. We will optimize the classic $k$-d tree
and use it to remove empty regions and increase the compressibility for each
AMR level (to be detailed in Section 3.3).

##  3 Our Proposed Design

In this section, we propose a pre-processing approach for AMR data to leverage
high-dimensional SZ lossy compression at each AMR level. Specifically, we
propose three pre-process strategies to mitigate the issue of irregular data
distribution. We further propose Shared Huffman Encoding (SHE) and integrate
it into the SZ compressor to further improve the compression performance for
AMR data. We also propose an adaptive approach to select the best-fit pre-
processing strategy based on the data density of each level.

###  3.1 Ghost-Shell Padding for High-density Data

To compress the AMR data in 3D, besides the aforementioned 3D baseline, we can
also compress each level separately in 3D. In that way, however, the data will
be split into multiple levels, and each level will have many empty regions and
an irregular data distribution, as shown in Figure 3. A naive solution to
handle the irregular 3D data is to fill the empty regions with zeros and pass
a large 3D block to the compressor. Although the padded zeros will increase
the size of data for compression, for high-density data such as z10’s coarse
level shown in Figure 3(b) (i.e., about 77% density), the size overhead will
be small.

(a) z10 fine level

(b) z10 coarse level

Figure 3: Visualization of data distributions of an example AMR data “z10”,
where z = redshift. Non-empty regions are shown in red.

However, these padded zeros can also greatly reduce the performance of
compression, especially for prediction-based lossy compression such as SZ,
because these zeros can significantly affect the prediction accuracy of SZ,
resulting in high compression errors on the boundaries, as shown in Figure
5(a). More specifically, as mentioned in Section 3.2, SZ uses each point’s
neighboring points’ values to predict its value. Thus, for those boundary
points which are adjacent to padded zeros, SZ will involve zero(s) in the
prediction, while the actual values of these empty regions are typically non-
zeros (saved in other AMR levels), which will seriously mislead the
prediction.

Figure 4: A 2D example of GSP approach. Non-empty blocks are in navy blue;
padded blocks are in light blue/red; padded blocks based on more than one non-
empty neighbor are in red.

(a) ZF (CR=156.7, PSNR=32.8dB)

(b) GSP (CR=161.3, PSNR=33.5dB)

Figure 5: Visual comparison (one slice) of compression errors of two
approaches using SZ based on Nyx’s “baryon density” field (i.e., z10’s coarse
level, 77% density). Brighter means higher compression error. The error bound
is the relative error bound of $6.7\times 10^{-3}$.

Input: Data, $x$, $y$

Output: Data after padding

1 for  _each unit block $b_{i}$_ do

2 if  _$b_{i}$ is empty  and $b_{i}$ has non-empty neighbor _ then

3 for  _each non-empty neighbor $n_{j}$_ do

4 pad slice = avg (first $y$ slices of $n_{j}$ next to $b_{i}$);

5 if  _overlap edge_ then

6 $pad={pad}/{2}$;

7

8 else if  _overlap corner_ then

9 $pad={pad}/{3}$;

10

11 else

12 continue;

13

14 end if

15 add an $x$-layers pad slice to $b_{i}$ next to $n_{j}$;

16

17 end for

18

19 end if

20

21 end for

return padded Data

Algorithm 1 Proposed Ghost Shell Padding Method

To eliminate the above issue of padding zeroes, we propose to use a ghost-
shell padding strategy (GSP) to diffuse neighboring values to a padding layer.
Figure 4 illustrates the high-level idea, and the detailed algorithm is
described in Algorithm 1. Specifically, we first partition the data into unit
blocks and then pad each empty unit block by using the average of its non-
empty neighbors’ boundary data values.

Note that some empty unit blocks have more than one non-empty neighbor such as
the red box shown in Figure 4. For these blocks, we will use the average value
of all its neighbors for padding. Correspondingly, we will remove these padded
values in the decompression based on the saved padding information. Note that
since the padding process is only for non-empty blocks, this metadata overhead
is almost negligible for high-density data (e.g., 0.1%).

After padding, each boundary point will be predicted using the average of all
the boundary data in the unit block(s) to which it belongs or is neighbored.
As shown in Figure 5, compared to the zero-filling (ZF) approach, GSP can
significantly reduce the overall compression error, especially for the
boundary data. Moreover, the GSP approach can provide a similar compression
ratio to the ZF approach on this high-density data and hence a better rate-
distortion. A detailed evaluation will be presented in Section 4.

###  3.2 Optimized Sparse Tensor Representation for Low-density Data

When most of the regions in the data are empty (e.g., about 77% of the data is
empty in Figure 3(a)), the large amount of padded data would greatly increase
the size of data for compression, resulting in a low compression ratio.

To solve this issue, we propose to use a naive sparse-tensor-based approach
(called NaST) to remove the empty regions, as shown in Figure 6. NaST includes
four main steps in the compression process: (1) partition the 3D data into
multiple unit blocks, (2) remove the empty blocks, (3) linearize the remaining
3D blocks into a 4D array, and (4) pass the 4D array to the compressor. Note
that in the decompression process, we will put the unit blocks from the
decompressed 4D array back into the original data.

Figure 6:  Workflow of the naive sparse tensor (NaST) method (empty regions
marked in pink and non-empty regions marked in blue).

However, in order to completely remove the empty regions to form a sparse
representation, the unit block size needs to be relatively small compared to
the input data size (e.g., $16^{3}$ vs. $512^{3}$), resulting in a high
proportion of data on the boundary. While the linearized unit blocks are
usually not adjacent in the original data, so the boundary data between them
are not smooth. Thus, it is harder for prediction-based compressors such as SZ
to predict the boundary data values. As a result, the NaST method without
optimizing the boundary data would have low compression performance.

To address the above problem, we propose an optimized sparse tensor
representation (called OpST) to effectively remove the empty regions as well
as maintain a relatively large unit block size so as to reduce the portion of
boundary data. A detailed description of our algorithm can be found in
Algorithm 2. We use a 2D example to demonstrate our approach, as illustrated
in Figure 7. Specifically, (1) we partition the data into many small unit
blocks. (2) For each unit block, we use the dynamic programming method to
initiate an array $BS$ to save the dimension/size of the maximum square whose
bottom-right corner is that unit block (line 6, which will be discussed in the
next paragraph). (3) We extract the sub-blocks (composed of multiple unit
blocks) from the original data according to the sizes saved in $BS$ (lines
13). (4) Since the original data will be changed after the extraction, we need
to partially update $BS$ based on maxSide (lines 14, will be discussed later).
We loop (3) and (4) from the bottom-right corner to the top-left corner until
the original data is empty. (5) After extracting all the sub-blocks, we put
them into multiple 3D arrays (to be compressed) based on their sizes. Note
that the sub-blocks with the same size will be merged into the same array for
easy compression.

Figure 7: A 2D example of our proposed OpST approach. The sub-blocks are
extracted according to our optimized sizes saved in $BS$. E.g., a 2-by-2 sub-
block $B_{0}$ is extracted according to $BS_{1}[2][1]$.

Input: Sparse 3D data S

Output: multiple 4D array $D_{n}$

1 for  _each unit block $b(x,y,z)$_ do

2 if  _$b(x,y,z)$ is non-empty_ then

3 if  _x is 0 or y is 0  or z is 0_ then

4 $BS(x,y,z)=1$

5 else

$BS(x,y,z)=\min(BS(x-1,y,z),~{}BS(x,y-1,z),~{}BS(x,y,z-1),~{}BS(x-1,y-1,z),~{}BS(x,y-1,z-1),~{}BS(x-1,y,z-1),~{}BS(x-1,y-1,z-1))+1$
;

/* BS(x,y,z) is the dimension size of the maximum cube whose bottom right rear
corner is the unit block with index (x,y,z) in the original data */

6 $maxSide=\max(maxSide,~{}BS(x,y,z))$

7 end if

8

9 end if

10

11 end for

12for  _each unit block $b(x,y,z)$_ do

13 if  _$BS(x,y,z)\geq 1$_ then

$size=BS(x,y,z)$ $D_{size}\leftarrow
S((x-size:x)*blkSize,~{}(y-size:y)*blkSize,~{}(z-size:z)*blkSize)$ ;

/* put the sub-block to the according to 4D array */

14 $b(x-size:x,~{}y-size:y,~{}z-size:z)\leftarrow empty$
$BS(x-size:x,~{}y-size:y,~{}z-size:z)=0$
$BS=updateBs(BS,~{}x,~{}y,~{}z,~{}maxSide)$

15 end if

16

17 end for

return $D_{n}$

Algorithm 2 Proposed Optimized Sparse Tensor Method

When initializing the $BS$ in step (2), we start with the $b^{\prime}[i][j]$
with $i=0$ or $j=0$ (i.e., on the top-left edge), where
$b^{\prime}[\cdot][\cdot]$ are the unit blocks: if $b^{\prime}[i][j]$ is
empty, we will set $BS[i][j]$ to 0 otherwise 1. For the remaining unit blocks,
if it is empty, $BS[i][j]$ will be 0; otherwise, $BS[i][j]$ will be set to 1
plus the minimum value among its three neighboring blocks (i.e., upper block,
left block, and upper-left block). In other words, we have
$BS[i][j]=1+\min(BS[i][j-1],BS[i-1][j],BS[i-1][j-1])$ for the 2D case. For
example, $BS_{1}[2][1]$ is 2 because all its upper-left neighbors are 1 (as
shown in Figure 7). However, both $BS_{1}[1][1]$ and $BS_{2}[1][2]$ can only
reach 1 because one of their neighbors is set to 0, having no chance to form a
sub-block with the size of 2.

Moreover, as mentioned in step (3), we need to update $BS$ after each
extraction. Specifically, for each sub-block we extract, we have to set its
corresponding values in $BS$ to zeros. For instance, as shown in Figure 7,
after we extract a 2-by-2 sub-block $B_{0}$ at $BS_{1}[2][1]$, we need to set
$BS_{2}[1][0]$, $BS_{2}[1][1]$, $BS_{2}[2][0]$, and $BS_{2}[2][1]$ to zeros.
In addition, we also need to recalculate a part of $BS$ (line 17 in Algorithm
2) because the extraction could influence other $BS$ values. For example, we
need to recalculate $BS_{2}[1][2]$ (marked in bold orange) after extracting
$B_{0}$. Note that this update is a partial update as the $BS$ values to be
updated will be bounded by maxSide which is the dimension size of the largest
cube in the dataset (line 7).

Similar to NaST, in decompression, we will put the sub-blocks back to
reconstruct the data based on the saved coordinates. Note that after our
optimization, each sub-block size will be relatively large (e.g., $96^{3}$ vs
the original data size of $512^{3}$), the overhead of saving the coordinates
of all the sub-blocks will be negligible (e.g., 0.1%).

(a) NaST(CR=245, PSNR=77.5dB)

(b) OpST(CR=248, PSNR=78.0dB)

Figure 8: Visual comparison (one slice) of compression errors of two
approaches using SZ based on Nyx’s “baryon density” field (i.e., z10’s fine
level, 23% density). Brighter means higher compression error. The error bound
is the relative error bound of $7.2\times 10^{-4}$.

Finally, we show a visual comparison of the compression quality between NaST
and OpST in Figure 8. Note that both use SZ with the same error bound.
Brighter means more errors. We can observe that compared to the NaST method,
OpST can significantly reduce the overall compression error, especially for
the data points on the boundary. It is worth noting that even with a lower
error, our OpST can still provide a higher compression ratio than NaST. This
is because our proposed optimization will generate larger sub-blocks, which
provide more information for prediction-based lossy compressors such as SZ to
achieve better rate-distortion. A detailed evaluation will be shown in Section
4.

###  3.3 Adaptive $k$-D Tree for Medium-density Data

The OpST approach proposed for low-density data, however, has a high
computation overhead, especially when the data is relatively dense. This is
because, on one hand, OpST needs to update BS based on maxSide for each
extraction of a sub-block, while the larger the maxSide, the more values in BS
that need to be updated; on the other hand, maxSide is the dimension size of
the largest non-empty cube in the dataset, which is highly related to the
density of the dataset. Thus, the time complexity of OpST can be expressed as
$O(N^{2}\cdot d)$, where $N$ is the unit block number and $d$ is the density.
Note that here density describes how dense the data is. For example, the
density of 77% means that 23% of the data is empty. Clearly, when the density
of an AMR level is relatively high, using OpST will be relatively time-
consuming.

Figure 9: 2D example of adaptive $k$-D tree. Sub-block will be adaptively
split to effectively remove empty regions and get bigger full sub-blocks.

Input: data block $d$, counts information

Output: $k$-d tree

1 node.count $\leftarrow$ counts information;

2 if  _$d$ is empty  or $d$ is full_ then

continue ;

/* stop splitting */

3

4else

5 if  _$d$ is a cube_ then

6 split d equally into 8 oct-blocks: $s_{1},\cdots,s_{8}$;

7 get the counts $c_{1},...c_{8}$ for $s_{1},\cdots,s_{8}$;

8 find the maxDiff partition $d_{1}$,$d_{2}$;

9 node.left = AKDTree ($d_{1}$, four $c_{i}$ of $d_{1}$);

10 node.right = AKDTree ($d_{2}$, four $c_{i}$ of $d_{2}$);

11

12 else if  _$d$ is a flat cuboid_ then

13 get the counts $c_{1},\cdots,c_{4}$ from counts information;

14 find the maxDiff partition $d_{1}$, $d_{2}$;

15 node.left = AKDTree ($d_{1}$, two $c_{i}$ of $d_{1}$);

16 node.right = AKDTree ($d_{2}$, two $c_{i}$ of $d_{2}$);

17

18 else if  _$d$ is a slim cuboid_ then

19 get the counts $c_{1},c_{2}$ from counts information;

20 split $d$ along the largest dimension to get $d_{1}$,$d_{2}$;

21 node.left = AKDTree ($d_{1}$, $c_{1}$);

22 node.right = AKDTree ($d_{2}$, $c_{2}$);

23

24

25 end if

return node;

Algorithm 3 Dynamic $k$-D Tree

To address the above high overhead issue of OpST, we propose an adaptive $k$-d
tree, called AKDTree, to remove empty regions and extract sub-blocks
(containing multiple unit blocks). AKDTree has a lower time complexity of
$O(\frac{1}{3}N\cdot\log N)$ (will be discussed later). Figure 9 shows a
simple 2D example. Specifically, (1) we partition the data into small unit
blocks. (2) We use a tree to hierarchically represent the whole data. Each
node in the tree is associated with a sub-block of the data. Moreover, each
node stores the number of non-empty unit blocks in the sub-block associated
with the node. (3) For each node, we split its associated sub-block from the
middle along one dimension to form two sub-blocks for its two children. Note
that we select one dimension which can maximize the difference of the numbers
of non-empty unit blocks of the two children (will be discussed in the next
two paragraphs). (4) We keep splitting a node until it has no empty unit block
or itself is empty. (5) Once finishing the construction of the tree, we
collect all the leaf nodes and send them to the compressor. Note that a non-
empty leaf node does not have any empty unit block; otherwise, it will keep
splitting. Thus, a leaf node must be an empty or full node, as shown in Figure
9. The detailed algorithm is described in Algorithm 3.

As mentioned in step (3), we are distributing the non-empty unit blocks
unevenly to two children for each node because we attempt to get as many leaf
nodes with large sub-block sizes as possible. If we keep splitting sub-blocks
in a fixed way, for instance, first split along the $x$-axis, second split
along the $y$-axis, third split along the $x$-axis, fourth split along the
$y$-axis, and so on, we will get a 2-by-2 sub-block for the node $n[2][2]$ as
shown in the dashed box, while its largest possible sub-block could be 4 by 2
as shown in Figure 9.

Figure 10: Example of adaptive splitting, different shapes will have a
different number of choices for splitting.

To select one of the dimensions to unevenly distribute its non-empty unit
blocks to the two children, we now present our dynamic splitting approach. We
categorize nodes into three different types: “cube” nodes, “flat” nodes, and
“slim” nodes, whose dimension ratios are 1:1:1, 2:2:1, 2:1:1, respectively.
First of all, for the cube node $d$, we first divide it into eight oct-blocks,
i.e., $s_{1}$, $s_{2}$, $\cdots$, $s_{8}$ (as shown in Figure 10), each sized
$\frac{n}{2}^{3}$. Here $n$ is the dimension size of the original data. Then,
we can get the counts of non-empty unit blocks of the eight oct-blocks, i.e.,
$c_{1}$, $c_{2}$, $\cdots$, $c_{8}$. After that, we will decide along which
dimension to split the cube node $d$ based on the counts. Specifically, we can
calculate the following three differences:

$\displaystyle\text{diff}_{x}$
$\displaystyle=|c_{1}+c_{3}+c_{5}+c_{7}-c_{2}-c_{4}-c_{6}-c_{8}|,$
$\displaystyle\text{diff}_{y}$
$\displaystyle=|c_{1}+c_{2}+c_{5}+c_{6}-c_{3}-c_{4}-c_{7}-c_{8}|,$
$\displaystyle\text{diff}_{z}$
$\displaystyle=|c_{1}+c_{2}+c_{3}+c_{4}-c_{5}-c_{6}-c_{7}-c_{8}|.$

Finally, we compare these three values and choose the dimension with the
maximum difference to split. For example, if the maximum difference is
$\text{diff}_{z}$, we will split $d$ along the z-axis (i.e., the pink 2D plane
shown in Figure 10) and get two flat nodes $d_{1}$ and $d_{2}$. For the flat
nodes such as $d_{1}$, we can reuse $c_{1}$, $\cdots$, $c_{4}$ to decide
whether to split $d_{1}$ along the x-axis or y-axis by choosing the larger one
among the following two differences.

$\displaystyle\text{diff}_{x}=|c_{1}+c_{3}-c_{2}-c_{4}|,\;\text{diff}_{y}=|c_{1}+c_{2}-c_{3}-c_{4}|.$

For the slim nodes such as $d_{11}$, we simply split it along the x-axis to
get two cube nodes $s_{1}$ and $s_{2}$. This process (i.e., cube
nodes$\to$flat nodes$\to$slim nodes) in step (3) will be looped until the node
becomes a leaf node (empty or full).

Note that based on the above description, the counting process is required for
every three nodes in each three path (i.e., only for the “cube” nodes). Thanks
to this dynamic splitting approach, we can lower the time complexity of the
AKDTree algorithm to $O\big{(}\frac{1}{3}\cdot N\cdot\log N\big{)}$, where $N$
is the number of unit blocks while extracting as many relatively large sub-
blocks without empty unit block as possible.

In addition, after the dynamic splitting, we will have a series of sub-blocks
with the same size but in different directions (e.g., 2:2:1, 2:1:2, 1:2:2). We
will align the sub-blocks with the same size based on their splitting
dimensions (instead of in-memory transposing them), merge them into an array,
and feed multiple merged arrays to the compression.

(a) d = 23

(b) d = 58

(c) d = 64

(d) d = 42

(e) d = 77

(f) d = 99

Figure 11: Compression performance comparison of GSP, OpST, and AKDTree on six
datasets with different densities.

###  3.4 Shared Huffman encoding

As mentioned in Section 1, 3.3 and 3.2, the OpST and AKDtree methods collect
and linearize the data blocks with the same shape into a 4D array and send it
to SZ. However, poor data locality/smoothness between non-neighbored data
blocks can negatively impact prediction accuracy, resulting in low data
quality. A potential solution could be compressing each data block
individually using SZ. However, this approach would result in low Huffman
encoding efficiency because the entire dataset would be divided into many
small blocks, requiring SZ to build a separate Huffman tree for each of these
small blocks. In other words, the original SZ method either requires
predicting and encoding small blocks together (by forcing them merged into 4D
arrays), leading to low prediction accuracy or predicting and encoding each
small block separately, which results in high Huffman encoding overhead.
Furthermore, even if the data blocks with the same shape can be compressed
together, the data blocks with different shapes still need to be compressed
separately, resulting in low Huffman encoding performance and a high time cost
of launching SZ multiple times.

Input: multiple data block $D_{1...n}$

Output: compressed data $S$

1 quantCode $Q$, regreCoeff $R$, compressed data $S$;

2

3for  _each data block $D_{i}$_ do

4 quantCode block $q_{i}$, regreCoeff block $r_{i}$;

5 $q_{i}$, $r_{i}$ = SZ.compress($D_{i}$);

6 $Q$.append($q_{i}$);

7 $R$.append($r_{i}$);

8

9 end for

10$S\leftarrow HuffmanEncode(Q)$;

$S\leftarrow HuffmanEncode(R)$;

/* end compression */

return S;

Algorithm 4 SZ compression with SHE

To this end, for OpST and AKDtree, we propose a shared Huffman encoding
technique to predict data blocks separately while encoding them together using
a single shared Huffman tree. The detail is shown in Algorithm 4. Each data
block is first predicted and quantized separately. Then, the quantization
codes and regression coefficients of each data block are aggregated to build a
shared Huffman tree and encoded at one time. This approach improves the
prediction performance of SZ without introducing high time overhead to the
encoding of SZ. As shown in Figure 3.4, compared to the original AKDtree,
AKDtree with SHE can significantly reduce the overall compression error,
especially for the data located at the boundary of data blocks, leading to
significant PSNR improvement, as shown in Figure 13. Furthermore, the use of
SHE reduces the number of Huffman trees needed for the data (since we do not
need separate Huffman trees for different block shapes), thereby improving
encoding efficiency and compression ratio.

(a) AKDtree
(CR=222, PSNR=78.5dB)

(b) AKDtree + SHE
(CR=231, PSNR=79.6dB)

Figure 12: Visual comparison (one slice) of compression errors of two
approaches using SZ based on Nyx’s “baryon density” field (i.e., z10’s fine
level, 23% density). Brighter means higher compression error. The error bound
is the relative error bound of $4.8\times 10^{-4}$. Figure 13: Comparison of
original AKDTree and AKDTree with SHE on Nyx’s “baryon density” field (i.e.,
z10’s fine level, 23% density).

###  3.5 Hybrid Compression Strategy

In this section, we propose a solution to adaptively choose a best-fit
compression strategy from our proposed OpST with SHE (OpST+), AKDTree with SHE
(AKDTree+), and GSP based on the data characteristics (i.e., data density).
According to Section 3.2 and 3.3, OpST+ is more suitable for sparse (i.e.,
low-density) data, while AKDTree+ is designed to address the high time
overhead of OpST+ when the density of data increases. Thus, there should be a
data-density threshold to determine when to use OpST+ or AKDTree+.

Figure 14: Time overhead comparison of OpST and AKDTree on different datasets
with different densities.

To decide the threshold $T$ for switching between OpST+ and AKDTree+, we
perform a series of experiments, as shown in Figure 11. The figure shows that
OpST+ and AKDTree+ have almost identical compression performance in terms of
bit-rate and PSNR on all six datasets/levels (from different timesteps) with
different densities. Moreover, Figure 14 shows the time costs of OpST+ and
AKDTree+ (excluding compression). The figure demonstrates that the time of
AKDTree+ is relatively stable, while the time of OpST+ increases linearly with
the increase of data density. Overall, the only criterion for selecting OpST
or AKDTree is the time cost rather than the compression performance. This is
consistent with our previous design aim, that is, AKDTree is mainly designed
to address the high time overhead issue of OpST+. Since OpST+ and AKDTree+
have a similar speed when the density is around 50%, we propose to choose $T$
= 50 for choosing OpST+ or AKDTree+.

According to Section 3.1, GSP is designed to effectively handle the AMR level
with high data density to prevent the negative impact of data partitioning
without the use of SHE. In contrast, the data partition methods such as
AKDTree require small blocks to be compressed together without SHE, which can
significantly decrease prediction accuracy, as shown in Section 3.4. However,
the negative impact on prediction accuracy caused by the partition can be
eliminated by using SHE to compress each small block produced by the
partition, which incurs little overhead, while GSP introduces significant
overhead to the data size. As shown in Figure 11, OpST+ and AKDTree+
outperform GSP across all the densities. As a result, the improved partition
strategies using SHE can be a viable alternative to GSP for all levels.

In summary, our proposed hybrid compression approach is described as follows.

  1. 1.

When the density is smaller than $T=50\%$, we use OpST+ to remove empty
regions and then compress;

  2. 2.

When the density is larger than $T=50\%$, we use AKDTree+ to remove empty
regions and then compress.

##  4 Experimental Evaluation

###  4.1 Experimental Setup

Test data. Our evaluation mainly focuses on the AMReX framework [3],
particularly the Nyx cosmology simulation [9]. Nyx is a state-of-the-art
extreme-scale cosmology code using AMReX, which generates six fields including
baryon density, dark matter density, temperature, and velocities ($x$, $y$,
and $z$). We use nine datasets generated by three real-world simulation runs
with different numbers of AMR levels, simulating a region of 64 megaparsecs
(Mpc). For this data, Z is equal to the redshift, i.e. the displacement
distant galaxies and celestial objects, as seen in Tab I.

Specifically, the first run has two levels of refinement, with the coarse
level of $256^{3}$ grids and the fine level of $512^{3}$ grids. We’ve
collected four timesteps with the finest level density from 23% to 64%. The
second run has a maximum of four levels of refinement. It was initially
configured at the resolution of $128^{3}$ and gradually refined to $1024^{3}$.
This run collected three timesteps with the finest-level resolution of
$256^{3}$ (two levels), $512^{3}$ (three levels), and $1024^{3}$ (four
levels), respectively. The density of the finest level varies from 0.2% to
0.003%. The third run has three levels of refinement with the grid sizes of
$128^{3}$ at the coarsest level, $256^{3}$ at the intermediate level, and
$512^{3}$ at the finest level. The density of the finest level ranges from
0.87% to 0.90%.

Note that the density of the finest level describes how much of the data in
the dataset is at the highest resolution; a higher density of the finest level
means that more data is refined to the highest resolution. Usually, the data
density is gradually increasing at the finest level, within a single run.

Evaluation platform and compressor. The test platform is equipped with two
28-core Intel Xeon Gold 6238R processors and 384 GB of memory. We use SZ3
[40], which is easy to (de)couple lossless encoding from SZ compression, due
to its high modularity.

TABLE I: Our tested datasets.

Dataset | # Levels |  | Grid Size of Each Level
---
(Fine to Coarse)
| Density of Each Level
---
(Fine to Coarse)
Nyx_Run1_Z10 | 2 | 512, 256 | 23%, 77%
Nyx_Run1_Z5 | 2 | 512, 256 | 58%, 42%
Nyx_Run1_Z2 | 2 | 512, 256 | 63%, 37%
Nyx_Run2_T2 | 2 | 256, 128 | 0.2%, 99.8%
Nyx_Run2_T3 | 3 | 512, 256, 128 | 0.02%, 0.56%, 99.42%
Nyx_Run2_T4 | 4 | 1024, 512, 256, 128 | 3E-5, 0.02%, 2.2%, 97.7%
Nyx_Run3_Z1.5 | 3 | 512, 256, 128 | 0.87%, 13.88%, 85.25%
Nyx_Run3_Z1 | 3 | 512, 256, 128 | 0.90%, 14.70%, 84.40%

Comparison baselines. As discussed in Section 2, we have three 1D or 3D
comparison baselines. Specifically, (1) the 1D baseline (naive): each AMR
level is compressed separately as a 1D array; (2) the 1D baseline (zMesh)
[27]: we refer readers to Section 2 for more details about how the zMesh
approach reorganize the AMR data for 1D compression; and (3) the 3D baseline:
Different AMR levels are unified to the same resolution for 3D compression;
(4) the TAC baseline: each AMR level is compressed using OpST (without SHE),
AKDTree (without SHE), and GSP based on the data density. For more
information, we refer readers to [41] for more details.

(a) Run1_Z10 (finest-level density = 23%)

(b) Run1_Z5 (finest-level density = 58%)

(c) Run1_Z2 (finest-level density = 63%)

Figure 15:  Rate-distortion comparison of our approach and baselines on the
early time-step (Z10) to the late time-step (Z2) from Nyx run1.

###  4.2 Evaluation Metrics

We will evaluate the compression performance based on the following metrics:
(1) compression ratio or bit-rate (generic), (2) distortion quality (generic),
(3) compression throughput (generic), (4) rate-distortion (generic), (5) power
spectrum (cosmology specific), (6) Halo finder (cosmology specific).

Metric 1: To evaluate the size reduction as a result of the compression, we
use the compression ratio, defined as the ratio of the original data size
compared to the compressed data size, or bit-rate (bits/value), representing
the amortized storage cost of each value. For single/double floating-point
data, the bit-rate is 32/64 bits per value before compression. The compression
ratio and bit-rate have a mathematical relationship as their multiplication is
32/64 so that a lower bit rate means a higher compression ratio.

Metric 2: Distortion is another important metric used to evaluate lossy
compression quality. We use the peak signal-to-noise ratio (PSNR) to measure
the distortion quality.

$\textstyle\text{PSNR}=20\cdot\log_{10}\left(R_{X}\right)-10\cdot\log_{10}\left(\sum_{i=1}^{N}{e_{i}^{2}}/{N}\right),$

where $e_{i}$ is the difference between the original and decompressed values
for the point $i$, $N$ is the number of points, and $R_{X}$ is the value range
of $X$. Higher PSNR less error.

Metric 3: (De)compression throughputs are critical to improving the I/O
performance. We calculate the throughput based on the original data size and
(de)compression time.

Metric 4: Similar to prior work [14, 15, 42, 32, 20, 43, 35], we plot the
rate-distortion curve to compare the distortion quality with the same bit-
rate, for a fair comparison between different compression approaches, taking
into account diverse compression algorithms.

Metric 5: Matter distribution in the Universe has evolved to form
astrophysical structures on different physical scales, from planets to larger
structures such as superclusters and galaxy filaments. The two-point
correlation function $\xi(r)$, which gives the excess probability of finding a
galaxy at a certain distance $r$ from another galaxy, statistically describes
the amount of the Universe at each physical scale. The Fourier transform of
$\xi(r)$ is called the matter power spectrum $P(k)$, where $k$ is the comoving
wavenumber. The matter power spectrum describes how much structure exists at
each physical scale. We run the power spectrum on the baryon density field by
using a cosmology analysis tool called Gimlet. We compare the power spectrum
$p^{\prime}(k)$ of decompressed data with the original $p(k)$ and accept a
maximum relative error within 1% for all k < 10.

Metric 6: Halo finder aims to find the halos (over-densities) in the dark
matter distribution and output the positions, the number of cells, and the
mass for each halo it finds, respectively. Specifically, the halo-finder
algorithm [44] searches for the halos from all the simulated data, with the
following two criteria: (1) the mass of a data point must be greater than a
threshold (e.g., 81.66$\times$ of the average mass of the whole dataset) to
become a halo cell candidate [20, 43, 45], and (2) there must be enough halo
cell candidates in a certain area to form a halo. For decompressed data, some
of the information (mass and cells of halos) can be distorted from the
original.

(a) Run2_T2 (finest-level density = 0.2%)

(b) Run2_T3 (finest-level density = 0.02%)

(c) Run2_T4 (finest-level density = 3E-5)

Figure 16: Rate-distortion comparison of TAC+ (top-left) and baselines on
different time-steps from run2.

###  4.3 Evaluation on Rate-distortion

We first evaluate the rate-distortion of our proposed compression approach and
compare it with the baselines on different datasets. As shown in Figure 15, 16
and 17, our new approach with SHE, TAC+ (represented by the pink curve) yields
better performance than the original TAC (the yellow curve) without SHE. Also,
for the 1D baseline, our approach (top-left curve) outperforms the 1D baseline
across all 8 datasets. Furthermore, the performance of our approach is more
stable (i.e., smoother curve) than the 1D baseline. We can also find that
zMesh is slightly worse than the 1D baseline on our tested data as shown in
Figure 15, which will be explained in the next section.

For the 3D baseline, we observe that our solution has much better performance
when the finest level has a relatively low density or the decompressed data
has a high PSNR, as shown in Figure 15, 16 and 17. However, our approach
cannot dominate the 3D baseline as shown in Figure 15(a), 15(b), and 16(a),
when the following criteria are satisfied: (1) the AMR data has only two
levels of refinement, (2) the finest level has a relatively low density, and
(3) the decompressed data has low PSNR/bitrate.

We also find that as the finest level density increases, the bitrate threshold
where the 3D baseline performs better than TAC+ is increasing, while the
advantage of the 3D baseline over TAC+ is decreasing, as shown in Figure 15
and 16(a). Specifically, in Figure 16(a) (the finest level density is 0.2%),
TAC+ outperforms the 1D baseline when the bit-rate is above 1.2; in Figure
15(a) (the finest level density is 23%), the intersection is the bit-rate of
1.25; as the finest level density continues to grow up to 58% in Figure 15(c),
TAC+ is slightly worse than the 3D baseline until the bit-rate is larger than
1.5; when the finest level density reaches 63% in Figure 15(c), TAC+ and the
3D baseline yield (almost) the same performance.

(a) Run3_Z1.5 (density for each level: 0.87%, 13.88%, 85.25% )

(b) Run3_Z1 (density for each level: 0.90%, 14.70%, 84.40%)

Figure 17:  Rate-distortion comparison of our approach and baselines on
different time-steps from Nyx run3.

In the next section, we will discuss how the number of AMR levels, the density
of the finest level, and the bit-rate/PSNR of decompressed data affect the
performance of the 3D baseline and TAC+ in detail.

###  4.4 Discussion on Comparison with Baselines

On compression, zMesh is meant to improve the smoothness of the patch-based
AMR datasets by taking advantage of the data redundancy between each AMR level
(as described in the introduction). Thus, zMesh cannot improve the smoothness
if there is no data redundancy in the tree-structured AMR datasets (i.e., our
tested datasets). A simple example is used to illustrate this in Figure 18(b),
where the finer-level data has higher values because a grid will be refined
only if its value is larger than a certain threshold. For block-based AMR,
when a grid needs to be refined because of its high value, the value will
still remain in the level, resulting in a redundant value saved (i.e., the red
8). If one uses the original z-ordering to traverse the data level-by-level
(shown in Figure 18(b)), the reordered data will have three significant value
changes (i.e., from 2 to 8, from 8 to 1, and from 1 to 9). To solve this
issue, zMesh traverses the two AMR levels together based on the layout of the
2D array. The reordered data are “1-2-8-9-8-7-8-1”, which only has two
significant value changes (i.e., from 2 to 8 and from 8 to 1). Thus, zMesh can
improve the smoothness of patch-based AMR data.

However, as shown in Figure 18(a), for tree-structured AMR data (without
saving a redundant “8”), compared to the 1D baseline that compresses each
level separately, zMesh introduces two significant data changes (i.e., from 2
to 9 and from 8 to 1) as it traverses between two AMR levels. This explains
why zMesh is slightly worse than the 1D baseline on our data.

(a) Tree-based AMR data

(b) Patch-based AMR data

Figure 18: An example of how the 1D baseline, zMesh, and original z-order
reorder a simple 2D AMR data without and with redundancy. Orange: coarse
level, blue: fine level, red: redundant data.

When considering the 3D baseline, we observe that it works slightly better
than TAC+ in the following circumstances: (1) the AMR data has only two levels
of refinement, (2) the decompressed data has a low PSNR/bit-rate, and (3) the
finest level of the data has a relatively low density. We also observe that
(4) when the finest level density increases, the bit-rate intersection point
where the 3D baseline performs better than TAC+ shifts to a higher value, (5)
while the advantage of the 3D baseline over TAC+ becomes smaller.

We now explain each observation point individually. As mentioned in Section
2.4, the main disadvantage of the 3D baseline is the redundant data generated
by the up-sampling process, which becomes more significant for the AMR data
with more than two levels of refinement, resulting in poor compression
performance (point 1). In contrast, the 3D baseline has better data locality
and smoothness over TAC+ because it compresses the entire dataset together
without partitioning the data. This leads to better compressibility for the 3D
baseline, resulting in faster bitrate reduction and slower PSNR degradation
when the error bound increases, compared to TAC+. As a result, the 3D baseline
performs better at a low bit-rate (point 2).

However, when the finest level density increases, the smoothness advantage of
the 3D baseline will decrease. This is because TAC+ will have more data
located in the finest level, and the partition overhead of TAC+ is lower at
the finest level compared to coarser levels due to the larger (unit) block
size. This explains why the 3D baseline only outperforms TAC+ for AMR datasets
with low density (points 3 & 5). Also, note that for the dataset with low
density in the finest level, the overhead of redundant up-sampled data will be
high for the 3D baseline, especially for high bit-rates (because the 3D
baseline only has better smoothness in low bit-rates). Thus, the compression
performance of the 3D baseline will decrease rapidly as we aim for relatively
high PSNRs/bit-rates, leading to a smaller intersection bit-rate between the
3D baseline and TAC+. Clearly, as the finest level density increases, the
overhead of the redundant data will be smaller, resulting in a larger
intersection bit-rate between the 3D baseline and TAC+ (point 4).

###  4.5 Post-analysis Quality with Adaptive Error Bound

When factoring level-wise compression, our approach can apply different error
bounds to different AMR levels based on (1) the post-analysis metrics, (2) the
up-sampling rates of coarse levels, and (3) the rate-distortion trade-off
between different AMR levels. We now evaluate our approach with the two
cosmology-specific post-analysis metrics (i.e., power spectrum and halo
finder) to demonstrate the benefit of the adaptive error bound method. We
choose the dataset run1-Z2 for evaluation because TAC+ has a similar
performance as the 3D baseline on this dataset.

TABLE II: Overall compression/decompression throughput (MB/s) of different
approaches with different absolute error bounds.

EB_abs | Run1_Z10 | Run1_Z2 | Run2_T2 | Run2_T4 | Run3_Z1
---|---|---|---|---|---
1D | 3D | TAC | TAC+ | 1D | 3D | TAC | TAC+ | 1D | 3D | TAC | TAC+ | 1D | 3D | TAC | TAC+ | 1D | 3D | TAC | TAC+
1E+08 | 51 | 37 | 79 | 77 | 51 | 79 | 85 | 84 | 71 | 15 | 71 | 66 | 49 | 0.32 | 24 | 24 | 49 | 5.6 | 63 | 63
1E+09 | 81 | 48 | 94 | 93 | 55 | 92 | 103 | 97 | 79 | 21 | 85 | 79 | 53 | 0.39 | 26 | 26 | 53 | 6.5 | 78 | 75
1E+10 | 89 | 53 | 107 | 102 | 58 | 100 | 111 | 104 | 86 | 23 | 98 | 87 | 84 | 0.41 | 28 | 27 | 57 | 7.1 | 91 | 84

Figure 19 shows the motivation of performing rate-distortion trade-off between
different AMR levels. As the error bounds for the fine and coarse levels
increase, their bit rates will converge to a similar value. This means that
when the error bound is relatively large, the reduction in data size will be
insignificant compared to the compression error increment (i.e., the slopes of
both curves are very small). Therefore, we can say that when the error is
large, it is not worth trading data quality for size reduction.

Figure 19: Bit-rates with different error bounds using SZ lossy compression
for fine and coarse levels on Run1_Z2 dataset.

(a) Power spectrum error under 1% limit

(b) Power spectrum error under 0.01% limit

Figure 20: Power spectrum error (in relative) of the 3D baseline and TAC+
(uniform error bound) and TAC+ (adaptive error bound) on baryon density field
on run1-Z2. The red and blue dashed line is the 1% and 0.01% limit of
acceptable power spectrum error.

Power Spectrum Figure 20(a) shows that under the (almost) same compression
ratio, TAC+ (with the uniform error bound) has a better power-spectrum error
compared to the 3D baseline. Also note that TAC+ yields nearly lossless power-
spectrum distortion (less than 0.01%) under the compression ratio of
42$\times$, as shown in Figure 20(a)’s gray curve.

Now, let us follow the three steps mentioned at the beginning of this section
to adjust the error bound for each AMR level. First, the post-analysis
metric–power spectrum—needs to be run on the uniform-resolution data and
focuses on the global quality of data. Thus, the ideal error-bound
configuration/ratio for the fine and coarse levels on the uniform-resolution
data would be 1:1.

As aforementioned, the coarse level of the AMR dataset needs to be up-sampled
to uniform the resolution. As a result, the compression error of the coarse
level will be up-sampled as well, resulting in more error in the post-
analysis. Thus, we then need to give the coarse level a smaller error bound
based on the up-sample rate. Here the up-sample rate for Z2’s coarse level is
$2^{3}$, leading to an ideal error-bound ratio of the fine and coarse levels
changed to 8:1.

Finally, this 8:1 ratio needs to be adjusted based on the rate-distortion
trade-off as aforementioned. As shown in Figure 19, when using the error-bound
ratio of 8:1 (e.g., 4E+9 for the fine level and 5E+8 for the coarse level),
the error bound of the fine level is too large, resulting in an ineffective
rate-distortion trade-off. Thus, we can balance two levels by increasing the
error bound for the coarse level (to gain compression ratio) and decreasing
the error bound for the fine level (to add compression error), which can
achieve an overall rate-distortion benefit. Based on our experiments, we
adjust the error-bound ratio from 8:1 to 3:1, As shown in Figure 20(b), TAC+
with adaptive error bound can significantly improve the compression ratio
compared with uniform error bound under similar power spectrum error.

Halo finer We evaluate the mass change, and the number of cells change for the
three largest halos identified using the 3D baseline, TAC+ (with uniform error
bound), and TAC+ (with adaptive error bound), as shown in Table III. We can
see that TAC+ with uniform error bound produces better halo-finer analysis
quality than the 3D baseline.

Similar to the error-bound configuration analysis done for the power spectrum,
let us now adjust the error-bound ratio between the fine and coarse levels for
halo finder. The halo-finer analysis also requires uniform-resolution data as
input. However, different from the power-spectrum analysis, the halo-finder
analysis focuses more on high-value points at the fine level, since only high-
value data points qualify as halo candidates, as described in Section 4.2.
Note that this does not mean we can directly discard the coarse-level data
with small values as they still contribute to the average value of the
dataset, which is also an important parameter for the halo finder [44].
Therefore, we set the ideal error-bound ratio to 1:2 (i.e., fine level vs
coarse level) for the uniform-resolution data based on our massive
experiments. After that, considering the up-sampling rate of $2^{3}$, the
error-bounded ratio is changed to 4:1. Finally, we adjust the ratio to 2:1
based on the rate-distortion trade-off. Table III shows that TAC+ with
adaptive error bound obtains the minimal differences of the mass and cell
numbers.

TABLE III: Halo finder analysis with different methods.

| CR | Avg Rel Mass Diff | Avg Rel Cells Diff
---|---|---|---
3D baseline | 188.7 | 2.7E-04 | 2.2E-03
TAC+ (1:1) | 189.1 | 2.2E-04 | 2.1E-03
TAC+ (2:1) | 192.5 | 1.1E-04 | 9.0E-04

###  4.6 Evaluation on Time Overhead

We evaluate the overall throughput (including preprocessing and
(de)compression) on the datasets with different error bounds. Compared to the
3D baseline, the throughput of TAC+ is up to $11.8\times$ higher on the Run3
dataset, $75\times$ higher on the Run2 datasets, and $2\times$ higher on the
Run1 datasets. This is because the Run2 and Run3 datasets have a lower density
than the Run1 datasets at the finest level, resulting in a higher overhead of
redundant data for the 3D baseline, which is consistent with our discussion in
Section 4.4.

We note that TAC+ performs better than the 1D baseline on all the tested
datasets except the T4 datasets, due to its relatively long data-partition
time compared to the total time on the small-sized datasets. TAC+ has almost
the same throughput as TAC. The very slight decrease is because TAC+
compresses the data in a more fine-grained manner. We exclude zMesh in this
evaluation as it is theoretically slower than the 1D baseline due to the extra
z-ordering and provides worse rate-distortion according to our evaluation.

##  5 Conclusion and Future Work

In conclusion, this paper leverages 3D compression for AMR data on a systemic
level. We propose three pre-processing strategies that can adapt based on the
density of each AMR level. Our approach improves the compression ratio
compared to the state-of-the-art approach by up to 4.9$\times$ under the same
data quality loss. With our level-wised compression approach, we are able to
tune the error-bound ratio of fine and coarse levels to be 3:1 and 2:1 for
better power-spectrum and halo-finder analyses, respectively, under the same
compression ratio. In future work, we will apply TAC+ to more AMR simulations
and improve its throughput on multi-core CPUs using OpenMP and GPUs using
CUDA.

## References

  * [1] C. Burstedde, O. Ghattas, G. Stadler, T. Tu, and L. C. Wilcox, “Towards adaptive mesh pde simulations on petascale computers,” _Proceedings of Teragrid_ , vol. 8, 2008.
  * [2] A. Dubey, A. Almgren, J. Bell, M. Berzins, S. Brandt, G. Bryan, P. Colella, D. Graves, M. Lijewski, F. Löffler _et al._ , “A survey of high level frameworks in block-structured adaptive mesh refinement packages,” _Journal of Parallel and Distributed Computing_ , vol. 74, no. 12, pp. 3217–3227, 2014.
  * [3] W. Zhang, A. Almgren, V. Beckner, J. Bell, J. Blaschke, C. Chan, M. Day, B. Friesen, K. Gott, D. Graves _et al._ , “Amrex: a framework for block-structured adaptive mesh refinement,” _Journal of Open Source Software_ , vol. 4, no. 37, pp. 1370–1370, 2019.
  * [4] J. M. Stone, K. Tomida, C. J. White, and K. G. Felker, “The athena++ adaptive mesh refinement framework: Design and magnetohydrodynamic solvers,” _The Astrophysical Journal Supplement Series_ , vol. 249, no. 1, p. 4, 2020\.
  * [5] A. S. Almgren, J. B. Bell, M. J. Lijewski, Z. Lukić, and E. Van Andel, “Nyx: A massively parallel amr code for computational cosmology,” _The Astrophysical Journal_ , vol. 765, no. 1, p. 39, 2013.
  * [6] B. Runnels, V. Agrawal, W. Zhang, and A. Almgren, “Massively parallel finite difference elasticity using block-structured adaptive mesh refinement with a geometric multigrid solver,” _Journal of Computational Physics_ , vol. 427, p. 110065, 2021.
  * [7] S. Whitman, J. Brasseur, and P. Hamlington, “Simulation of bluff-body stabilized flames with pelec, an exascale combustion code,” 2018.
  * [8] K. Sverdrup, N. Nikiforakis, and A. Almgren, “Highly parallelisable simulations of time-dependent viscoplastic fluid flow with structured adaptive mesh refinement,” _Physics of Fluids_ , vol. 30, no. 9, p. 093102, 2018.
  * [9] Nyx, https://github.com/AMReX-Astro/Nyx, 2021.
  * [10] P. Deutsch, “Gzip file format specification version 4.3,” 1996.
  * [11] Zstandard, http://facebook.github.io/zstd/, 2020.
  * [12] S. W. Son, Z. Chen, W. Hendrix, A. Agrawal, W.-k. Liao, and A. Choudhary, “Data compression for the exascale computing era-survey,” _Supercomputing Frontiers and Innovations_ , vol. 1, no. 2, pp. 76–88, 2014\.
  * [13] S. Di and F. Cappello, “Fast error-bounded lossy hpc data compression with sz,” in _2016 ieee international parallel and distributed processing symposium (ipdps)_. IEEE, 2016, pp. 730–739.
  * [14] D. Tao, S. Di, Z. Chen, and F. Cappello, “Significantly improving lossy compression for scientific data sets based on multidimensional prediction and error-controlled quantization,” in _2017 IEEE International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, 2017, pp. 1129–1139.
  * [15] X. Liang, S. Di, D. Tao, S. Li, S. Li, H. Guo, Z. Chen, and F. Cappello, “Error-controlled lossy compression optimized for high compression ratios of scientific datasets,” in _2018 IEEE International Conference on Big Data_. IEEE, 2018.
  * [16] P. Lindstrom, “Fixed-rate compressed floating-point arrays,” _IEEE Transactions on Visualization and Computer Graphics_ , vol. 20, no. 12, pp. 2674–2683, 2014.
  * [17] M. Ainsworth, O. Tugluk, B. Whitney, and S. Klasky, “Mgard: A multilevel technique for compression of floating-point data,” in _DRBSD-2 Workshop at Supercomputing_ , 2017.
  * [18] R. Ballester-Ripoll, P. Lindstrom, and R. Pajarola, “Tthresh: Tensor compression for multidimensional visual data,” _IEEE Transactions on Visualization and Computer Graphics_ , vol. 26, no. 9, pp. 2891–2903, 2020.
  * [19] F. Cappello, S. Di, S. Li, X. Liang, A. M. Gok, D. Tao, C. H. Yoon, X.-C. Wu, Y. Alexeev, and F. T. Chong, “Use cases of lossy compression for floating-point data in scientific data sets,” _The International Journal of High Performance Computing Applications_ , 2019.
  * [20] S. Jin, P. Grosset, C. M. Biwer, J. Pulido, J. Tian, D. Tao, and J. Ahrens, “Understanding gpu-based lossy compression for extreme-scale cosmological simulations,” in _2020 IEEE International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, 2020, pp. 105–115.
  * [21] P. Grosset, C. M. Biwer, J. Pulido, A. T. Mohan, A. Biswas, J. Patchett, T. L. Turton, D. H. Rogers, D. Livescu, and J. Ahrens, “Foresight: analysis that matters for data reduction,” in _SC20: International Conference for High Performance Computing, Networking, Storage and Analysis_. IEEE, 2020, pp. 1–15.
  * [22] T. Lu, Q. Liu, X. He, H. Luo, E. Suchyta, J. Choi, N. Podhorszki, S. Klasky, M. Wolf, T. Liu _et al._ , “Understanding and modeling lossy compression schemes on hpc scientific data,” in _2018 IEEE International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, 2018, pp. 348–357.
  * [23] A. H. Baker, H. Xu, J. M. Dennis, M. N. Levy, D. Nychka, S. A. Mickelson, J. Edwards, M. Vertenstein, and A. Wegener, “A methodology for evaluating the impact of data compression on climate simulation data,” in _Proceedings of the 23rd International Symposium on High-Performance Parallel and Distributed Computing_. ACM, 2014, pp. 203–214.
  * [24] A. H. Baker, H. Xu, D. M. Hammerling, S. Li, and J. P. Clyne, “Toward a multi-method approach: Lossy data compression for climate simulation data,” in _International Conference on High Performance Computing_. Springer, 2017, pp. 30–42.
  * [25] A. M. Gok, S. Di, Y. Alexeev, D. Tao, V. Mironov, X. Liang, and F. Cappello, “Pastri: Error-bounded lossy compression for two-electron integrals in quantum chemistry,” in _2018 IEEE international conference on cluster computing (CLUSTER)_. IEEE, 2018, pp. 1–11.
  * [26] X.-C. Wu, S. Di, E. M. Dasgupta, F. Cappello, H. Finkel, Y. Alexeev, and F. T. Chong, “Full-state quantum circuit simulation by using data compression,” in _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ , 2019, pp. 1–24.
  * [27] H. Luo, J. Wang, Q. Liu, J. Chen, S. Klasky, and N. Podhorszki, “zmesh: Exploring application characteristics to improve lossy compression ratio for adaptive mesh refinement,” in _2021 IEEE International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, 2021, pp. 402–411.
  * [28] G. K. Wallace, “The JPEG still picture compression standard,” _IEEE Transactions on Consumer Electronics_ , vol. 38, no. 1, pp. xviii–xxxiv, 1992\.
  * [29] D. Le Gall, “Mpeg: A video compression standard for multimedia applications,” _Communications of the ACM_ , vol. 34, no. 4, pp. 46–58, 1991.
  * [30] “2021 r&d 100 award winners,” https://www.rdworldonline.com/rd-100-2021-winner/sz-a-lossy-compression-framework-for-scientific-data/, 2021\.
  * [31] K. Zhao, S. Di, M. Dmitriev, T.-L. D. Tonellot, Z. Chen, and F. Cappello, “Optimizing error-bounded lossy compression for scientific data by dynamic spline interpolation,” in _2021 IEEE 37th International Conference on Data Engineering_. IEEE, 2021, pp. 1643–1654.
  * [32] X. Liang, Q. Gong, J. Chen, B. Whitney, L. Wan, Q. Liu, D. Pugmire, R. Archibald, N. Podhorszki, and S. Klasky, “Error-controlled, progressive, and adaptable retrieval of scientific data with multilevel decomposition,” in _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ , 2021, pp. 1–13.
  * [33] F. Wang, N. Marshak, W. Usher, C. Burstedde, A. Knoll, T. Heister, and C. Johnson, “Cpu ray tracing of tree-based adaptive mesh refinement data,” _Computer Graphics Forum_ , vol. 39, pp. 1–12, 06 2020.
  * [34] G. Harel, J.-B. Lekien, and P. Pébaÿ, “Two new contributions to the visualization of amr grids: I. interactive rendering of extreme-scale 2-dimensional grids ii. novel selection filters in arbitrary dimension,” 03 2017\.
  * [35] K. Zhao, S. Di, X. Liang, S. Li, D. Tao, Z. Chen, and F. Cappello, “Significantly improving lossy compression for hpc datasets with second-order prediction and parameter optimization,” in _Proceedings of the 29th International Symposium on High-Performance Parallel and Distributed Computing_ , 2020, pp. 89–100.
  * [36] J. Bentley, “Multidimensional binary search trees used for associative searching,” _communications of the ACM September, 1975. vol. 18: pp. 509-517 : ill. includes bibliography._ , vol. 18, 01 1975.
  * [37] D. Hoang, H. Bhatia, P. Lindstrom, and V. Pascucci, “High-quality and low-memory-footprint progressive decoding of large-scale particle data,” in _2021 IEEE 11th Symposium on Large Data Analysis and Visualization (LDAV)_. IEEE Computer Society, 2021, pp. 32–42.
  * [38] G. Cirio, G. Lavoué, and F. Dupont, “A framework for data-driven progressive mesh compression.” in _GRAPP_ , 2010, pp. 5–12.
  * [39] O. Devillers and P.-M. Gandoin, “Geometric compression for interactive transmission,” _Proc. Visualization ’00_ , 01 2000.
  * [40] X. Liang, K. Zhao, S. Di, S. Li, R. Underwood, A. M. Gok, J. Tian, J. Deng, J. C. Calhoun, D. Tao, Z. Chen, and F. Cappello, “Sz3: A modular framework for composing prediction-based error-bounded lossy compressors,” _IEEE Transactions on Big Data_ , pp. 1–14, 2022.
  * [41] D. Wang, J. Pulido, P. Grosset, S. Jin, J. Tian, J. Ahrens, and D. Tao, “Tac: Optimizing error-bounded lossy compression for three-dimensional adaptive mesh refinement simulations,” in _Proceedings of the 31st International Symposium on High-Performance Parallel and Distributed Computing_ , 2022, p. 135–147.
  * [42] X. Liang, S. Di, D. Tao, Z. Chen, and F. Cappello, “An efficient transformation scheme for lossy data compression with point-wise relative error bound,” in _2018 IEEE International Conference on Cluster Computing (CLUSTER)_. IEEE, 2018, pp. 179–189.
  * [43] S. Jin, J. Pulido, P. Grosset, J. Tian, D. Tao, and J. Ahrens, “Adaptive configuration of in situ lossy compression for cosmology simulations via fine-grained rate-quality modeling,” in _Proceedings of the 30th International Symposium on High-Performance Parallel and Distributed Computing_ , 2021, pp. 45–56.
  * [44] M. Davis, G. Efstathiou, C. S. Frenk, and S. D. White, “The evolution of large-scale structure in a universe dominated by cold dark matter,” _The Astrophysical Journal_ , vol. 292, pp. 371–394, 1985.
  * [45] B. Fang, D. Wang, S. Jin, Q. Koziol, Z. Zhang, Q. Guan, S. Byna, S. Krishnamoorthy, and D. Tao, “Characterizing impacts of storage faults on hpc applications: A methodology and insights,” in _2021 IEEE International Conference on Cluster Computing (CLUSTER)_. IEEE, 2021, pp. 409–420.

|  Daoce Wang is a PhD student in Intelligent Systems Engineering at Indiana
University Bloomington. He received his bachelor’s degree in Computer Science
from University of Electronic Science and Technology of China in 2018 and his
master’s degree in Computer Science from University of Florida in 2020. His
research interests include scientific data reduction, AMR simulations, and
parallel algorithms.
---|---
|  Jesus Pulido is a research scientist in the Data Science at Scale Team at
Los Alamos National Laboratory. Pulido received his Ph.D. in Computer Science
at University of California, Davis in 2019. Pulido specializes in data
analysis, data reduction, visualization, high performance computing, wavelets
and multi-resolution methods. He has experience in applications of image
sensors, astronomy, turbulence and cosmology.
---|---
|  Pascal Grosset is a scientist in the Data Science at Scale team at Los
Alamos National Laboratory. His primary research interests are large-scale
data analysis and visualization, and data reduction. He received his Ph.D. in
Computing: Graphics and Visualization from the University of Utah in 2016
where his research focused on large scale visualization.
---|---
|  Sian Jin is currently working toward the PhD degree majoring in computer
engineering at Indiana University Bloomington. His research interests include
High Performance Computing, Compression Algorithms, Artificial Neural
Networks, and Parallel Computing. He has published several papers in major
journals and international conferences including the SC, PPoPP, VLDB, HPDC,
IPDPS, and ICDE.
---|---
|  Jiannan Tian is a PhD candidate in Intelligent Systems Engineering at
Indiana University Bloomington. His research interests include lossy
compression for scientific data and error analysis, and GPU-centric computing.
His ongoing project including developing GPU-accelerated compression algorithm
and system design optimization of lossy compression.
---|---
|  Kai Zhao is an assistant professor of computer science at University of
Alabama at Birmingham. He received his Ph.D. in computer science from
University of California, Riverside in 2022. He received his bachelor’s degree
from Peking University in 2014. His research interests include high-
performance computing, scientific data management and reduction, and resilient
machine learning.
---|---
|  James P. Ahrens received the BS degree in computer science from the
University of Massachusetts at Amherst, Amherst, Massachusetts, in 1989, and
the PhD degree in computer science from the University of Washington, Seattle,
Washington, in 1996. He is a senior scientist with Applied Computer Science
Group, Los Alamos National Laboratory. His primary research interests include
visualization, computer graphics, data science, and parallel systems. He is
author of more than 100 peer reviewed papers and the founder/design lead of
ParaView, an open-source visualization tool designed to handle extremely large
data. ParaView is broadly used for scientific visualization and is in use at
supercomputing and scientific centers worldwide. He is the chair of the IEEE
Computer Society Technical Committee on Visualization and Graphics (VGTC).
---|---
|  Dingwen Tao is an associate professor at Indiana University Bloomington,
where he directs the High-Performance Data Analytics and Computing Lab. He
received his Ph.D. in Computer Science from University of California,
Riverside in 2018 and B.S. in Mathematics from University of Science and
Technology of China in 2013. He is the recipient of various awards including
NSF CAREER Award (2023), Amazon Research Award (2022), Meta Research Award
(2022), R&D100 Awards Winner (2021), IEEE Computer Society TCHPC Early Career
Researchers Award for Excellence in HPC (2020), NSF CRII Award (2020), and
IEEE CLUSTER Best Paper Award (2018). He is serving on the Technical Review
Board of IEEE Transactions on Parallel and Distributed Systems. He served as
the Program Co-chair of 2021 IEEE International Conference on Scalable
Computing and Communications and International Workshops on Big Data
Reduction. He is also a reviewer, program committee member, or session chair
of major HPC venues, such as SC, HPDC, ICS, IPDPS, CLUSTER, ICPP, CCGrid, and
HiPC. Email<EMAIL_ADDRESS>
---|---