# DuAT: Dual-Aggregation Transformer Network for Medical Image Segmentation

Feilong Tang  School of AI and Advanced Computing, Xi’an Jiaotong-Liverpool
University, Suzhou, China  Qiming Huang  School of AI and Advanced Computing,
Xi’an Jiaotong-Liverpool University, Suzhou, China  Jinfeng Wang  School of AI
and Advanced Computing, Xi’an Jiaotong-Liverpool University, Suzhou, China
Xianxu Hou  School of Computer Science and Software Engineering, Shenzhen
University, Shenzhen, China  Jionglong Su🖂 School of AI and Advanced
Computing, Xi’an Jiaotong-Liverpool University, Suzhou, China  Jingxin Liu🖂
School of AI and Advanced Computing, Xi’an Jiaotong-Liverpool University,
Suzhou, China

###### Abstract

Transformer-based models have been widely demonstrated to be successful in
computer vision tasks by modelling long-range dependencies and capturing
global representations. However, they are often dominated by features of large
patterns leading to the loss of local details (e.g., boundaries and small
objects), which are critical in medical image segmentation. To alleviate this
problem, we propose a Dual-Aggregation Transformer Network called DuAT, which
is characterized by two innovative designs, namely, the Global-to-Local
Spatial Aggregation (GLSA) and Selective Boundary Aggregation (SBA) modules.
The GLSA has the ability to aggregate and represent both global and local
spatial features, which are beneficial for locating large and small objects,
respectively. The SBA module is used to aggregate the boundary characteristic
from low-level features and semantic information from high-level features for
better preserving boundary details and locating the re-calibration objects.
Extensive experiments in six benchmark datasets demonstrate that our proposed
model outperforms state-of-the-art methods in the segmentation of skin lesion
images, and polyps in colonoscopy images. In addition, our approach is more
robust than existing methods in various challenging situations such as small
object segmentation and ambiguous object boundaries.

##  1 Introduction

Medical image segmentation is a computer-aided automatic procedure for
extracting the region of interest (RoI), e.g., tissues, lesions, and body
organs. It can assist clinicians by making diagnostic and treatment processes
more efficient and precise. For instance, colonoscopy is the gold standard for
detecting colorectal lesions, and accurately locating early polyps is of great
significance for clinical prevention of rectal cancer [17]. Likewise, melanoma
skin cancer is one of the most rapidly increasing cancers worldwide.
Segmentation of skin lesions from dermoscopic images is a critical step in
skin cancer diagnosis and treatment planning [33]. However, it is impractical
to manually annotate these structures in clinical practice due to the tedious,
time-consuming, and error-prone process. There is a growing need for automatic
as well as accurate image segmentation.

Figure 1:  Main challenges of medical images. Challenge1: Small object
segmentation are difficult to segment due to their low contrast and strong
camouflage. Challenge2: The object boundary for medical image is ambiguous due
to image acquisition influence.

With the rapid development of deep learning, an increasing number of deep
convolutional neural networks (DCNNs) [16, 18, 36, 54, 59] are proposed for
medical image segmentation. The limitations of the receptive field in DCNNs
make it difficult to capture the global representation. To solve this problem,
dilated convolution [11, 51, 53] is proposed for semantic segmentation task.
Furthermore, attention models [47, 21, 52] are developed to better capture
long-range context information. These alternatives achieve promising results
in semantic segmentation.

Transformer-based methods [1, 8, 9, 19, 40, 58, 46] have been proposed and
achieved comparable performance to state-of-the-art results. Transformer is
originally used for sequence-to-sequence predictive modelling in natural
language processing (NLP) tasks with pure attention structure that is good at
capturing long-range dependencies [42]. Vision Transformer-based methods [1,
8, 9, 19, 40, 58, 46] have also been proposed and demonstrate promising
performance. ViT [15] constructs a vector sequence by dividing each image into
fixed-size patches, subsequently applies a multi-head self-attention (MHSA)
and the Multilayer Perceptron (MLP) structure which demonstrate the advanced
learning ability for long-distance feature dependence. The recent works [14,
50, 45] demonstrate that the pyramid structure is applicable in Transformers
and more suitable for various downstream tasks. Unfortunately, capturing long-
range dependencies destroys part of local features, which could result in
overly smooth predictions for small objects and blurred boundaries between
objects.

Therefore, building a model that retains both local and global features
remains challenging. Li et al. [26] explore the local context for the
aggregated long-range relationship to be more accurately distributed in local
regions. Recently, some hybrid architectures of Transformer and CNN are
proposed [10, 35, 58], which aim to combine the advantages of both models.
Wang et al. [45] propose progressive locality decoder to emphasize local
features and restrict attention dispersion. Chen et al. [10] propose
TransUnet, which utilizes the underlying features of CNNs and subsequently
uses the transformers to model global interactions to strengthen local
features. However, feeding local information directly into the transformer
cannot precisely handle local context relationships, resulting in the local
information being overwhelmed by the dominant global context. Ultimately it
leads to inferior results in the medical image segmentation of small objects.

Considering the problem of unclear object boundary information in semantic
segmentation. Previous studies [27, 12, 14] explore fusing low-scale boundary
information and high-scale semantic information to better preserve boundary
details. Moreover, Zhang et al. [56] incorporate semantic information into
low-level features while embedding more spatial information into high-level
features to make up for the semantic and resolution gap between feature maps.
Takikawa et al. [39] and Zhen et al. [57] design a boundary stream and couple
the task of boundary and semantics modeling. Inspired by [56], we exploit the
boundary information to selectively fuse it with the high-level semantic
features, which will further enhance the semantic representations rather than
simply combining them.

In this paper, we propose a novel pyramid transformer for medical image
segmentation, referred to as the Dual-Aggregation Transformer Network (DuAT),
consisting of the Global-to-Local Spatial Aggregation (GLSA) to combine local
and global features, as well as the Selective Boundary Aggregation (SBA)
module to enhance the boundary information and locating the re-calibration
objects. We believe that global spatial features help locate large objects,
and local spatial features are crucial for identifying small ones. Finally,
boundary information is aggregated to fine-tune object boundaries and re-
calibrate coarse predictions, Specifically, the GLSA module is responsible for
extracting and fusing local and global spatial features from the backbone. We
separate the channels, one for global representation extracted by Global
context (GC) block [7], and the other for local information extracted by
multiple depth-wise convolutions. The SBA module aims to simulate the
biological visual perception process, distinguishing objects from background.
Specifically, it incorporates shallow- and deep-level features to establish
the relationship between body areas and boundary, enhancing the boundaries
characteristics. Our experimental results demonstrate three advantages of our
model: more lightweight, better learning ability and improved generalization
capability than previous state-of-the-art approaches.

In summary, the main contribution of this paper is three-fold.

  * $\bullet$

We propose a novel framework, named Dual-Aggregation Transformer Network
(DuAT), for medical images segmentation, which adapts the pyramid vision
transformer as encoder to extract more robust features than the existing CNN-
based methods.

  * $\bullet$

We design dual aggregation modules, Global-to-Local Spatial Aggregation (GLSA)
module and Selective Boundary Aggregation (SBA) module. Their purpose is to
solve two challenges that are intuitively demonstrated in Figure 1.
Specifically, the GLSA module simultaneously extract local spatial detail
information and global spatial semantic information, which reduces incorrect
information in the high-level features. The SBA module aims to fine-tune
object boundaries, which can well address the “ambiguous” problem of
boundaries.

  * $\bullet$

Extensive experiments on five polyp datasets (ETIS [43], CVC-ClinicDB [37],
CVC-ColonDB [3], EndoScene-CVC300 [38], Kvasir [22]), skin lesion dataset
(ISIC-2018 [13]) and 2018 Data Science Bowl [6] demonstrate that the proposed
DuAT methods advances the state-of-the-art (SOTA) performance.

Figure 2: The overall architecture of Dual-Aggregation Transformer Network
(DuAT). The entire model is divided into three parts: (a) pyramid vision
transformer (PVT) as backbone; (b) pyramid Global-to-Local Spatial Aggregation
(GLSA) Module; (c) Selective Boundary Aggregation (SBA) module and it shown on
the red box.

##  2 Related Work

###  2.1 Vision Transformer

Transformer has dominated the field of NLP with its MHSA layer to capture the
pure attention structure of long-range dependencies. Different from
convolutional layer, MHSA layer has dynamic weight and global receptive field,
which makes it more flexible and effective. Dosovitskiy et al. propose a
vision transformer (ViT) [15], which is an end-to-end model using the
Transformer structure for image recognition task. Specifically, it divides an
image into fixed-size patches, which are sequentially fed to multiple
Transformer encoder blocks to model the patches. In addition, previous work
has proved that the pyramid structure in convolutional networks is also
suitable for Transformer and various downstream tasks, such as Swin
Transformer [28], PVT [46], Segformer [50], etc. PVT requires less computation
than ViT and adopts the classical Semantic-FPN to deploy the task of semantic
segmentation.

In medical image segmentation, TransUNet [10] demonstrates that Transformer
can be used as powerful encoders for medical image segmentation. TransFuse
[55] is proposed to improve efficiency for global context modeling by fusing
transformers and CNNs. Furthermore, to train the model effectively on medical
images, Polyp-PVT [14] introduces Similarity Aggregation Module based on graph
convolution domain [31]. Inspired by these approaches, we propose a new
transformer-based medical segmentation framework, which can accurately locate
small objects.

Figure 3:  Overview of the Global-to-Local Spatial Aggregation Module GLSA, it
is composed of global spatial attention (GSA) and local spatial attention
(LSA).

###  2.2 Image Boundary Segmentation

Recently, learning additional boundary information has shown superior
performance in many image segmentation tasks. In the early research on FCN-
based semantic segmentation, Bertasius et al. [4] and Chen et al. [11] use
boundaries for post-starting to refine the result at the end of the network.
Recently, several approaches explicitly model boundary detection as an
independent sub-task in parallel with semantic segmentation for sharper
result. Ma et al. [32] explicitly exploit the boundary information for context
aggregation to further enhance the semantic representation of the model. Li et
al. [25] point out that the object boundary and body parts correspond to the
high-frequency and low-frequency information of an image, respectively, based
on which they decouple the body and edge with diverse supervisions. Ji et al.
[24] fuse the low-level edge-aware features and constraint it with the
explicit edge supervision. Different from the above works, we propose a novel
aggregation method to achieve more accurate localisation and boundary
delineation of objects.

##  3 Method

As given in Figure 2, our DuAT model consists of the following:- a pyramid
vision transformer (PVT) encoder, a SBA module, and GLSA module.

###  3.1 Transformer Encoder

Some recent studies [5, 50] report that vision transformers [15, 46] have
stronger performance and robustness to input disturbances (e.g., noise) than
CNNs. Inspired by this, we use the Transformer based on pyramid structure as
the encoder. Specifically, the pyramid vision transformer (PVT) [46] is
utilized as the encoder module for multi-level feature maps
{${\mathcal{F}_{i}|i\in(1,2,3,4)}$} extraction. Among these feature maps,
$\mathcal{F}_{1}$ gives detailed boundary information of target, and
$\mathcal{F}_{2}$, $\mathcal{F}_{3}$ and $\mathcal{F}_{4}$ provide high-level
features.

###  3.2 Selective Boundary Aggregation

As observed in [56, 25], shallow- and deep-layer features complement each
other. The shallow layer has less semantics but is rich in details, with more
distinct boundaries and less distortion. Furthermore, the deep level contains
a rich semantic information. Therefore, directly fusing low-level features
with high-level ones may result in redundancy and inconsistency. To address
this, we propose the SBA module, which selectively aggregate the boundary
information and semantic information to depict more fine-grained contour of
objects and the location of re-calibrate objects.

Different from previous fusion methods, we design a novel Re-calibration
attention unit (RAU) block that adaptively picks up mutual representations
from tow inputs ($F^{s},F^{b}$) before fusion. As given in Figure 2, the
shallow - and deep-level information is fed into the two RAU blocks by
different ways to make up for the missing spatial boundary information of the
high-level semantic features and the missing semantic information of low-level
features. Finally, the outputs of two RAU blocks are concatenated after a
$3\times 3$ convolution. This aggregation strategy realizes the robust
combination of different features and refines the rough features. The RAU
block function $PAU(\cdot,\cdot)$ process can be expressed as:

$\displaystyle
T_{1}^{\prime}=W_{\theta}(T_{1}),T_{2}^{\prime}=W_{\phi}(T_{2})$ (1)
$\displaystyle PAU(T_{1},T_{2})=T_{1}^{\prime}\odot T_{1}+T_{2}^{\prime}\odot
T_{2}\odot(\circleddash(T_{1}^{\prime}))+T_{1},$ (2)

where $T_{1},T_{2}$ are the input features, two linear mapping and sigmoid
functions $W_{\theta}(\cdot),W_{\phi}(\cdot)$ are applied to the input
features to reduce the channel dimension to 32 and obtain feature maps
$T_{1}^{\prime}$ and $T_{2}^{\prime}$. $\odot$ is Point-wise multiplication.
$\circleddash(\cdot)$ is the reverse operation by subtracting the feature
$T_{1}^{\prime}$, refining the imprecise and coarse estimation into an
accurate and complete prediction map [16]. We take a convolutional operation
with a kernel size of 1 $\times$ 1 as the linear mapping process. As a result,
the process of SBA can be formulated as:

$\displaystyle Z=C_{3\times 3}({\rm
Concat}(PAU(F^{s},F^{b}),PAU(F^{b},F^{s}))),$ (3)

where $C_{3\times 3}(\cdot)$ is a $3\times 3$ convolution with a batch
normalization and a ReLU activation layer.
$F^{s}\in\mathbb{R}^{\frac{H}{8}\times\frac{W}{8}\times 32}$ contains deep-
level semantic information after fusing the third and fourth layers from the
encoder, $F^{b}\in\mathbb{R}^{\frac{H}{4}\times\frac{W}{4}\times 32}$ is the
first layer with rich boundary details from the backbone. $Concat(\cdot)$ is
the concatenation operation along the channel dimension.
$Z\in\mathbb{R}^{\frac{H}{4}\times\frac{W}{4}\times 32}$ is the output of the
SBA module.

Figure 4: Qualitative results of different methods. (a) Inputs images, (b) GT,
which stands for the ground truths, (h) semantic segmentation maps produced by
our method, (c) U-Net [36], (d) U-Net++ [59], (e) PraNet [16], (f) TransUnet
[10], (g) Polyp-PVT [14].

###  3.3 Global-to-Local Spatial Aggregation

The attention mechanism strengthens the information related to the
optimization goal and suppresses irrelevant information. In order to capture
both global and local spatial features, we propose the GLSA module, which
fuses the results of two separate local and global attention units. As
demonstrated in Figure 3, this dual-stream design effectively preserves both
local and non-local modeling capabilities. Moreover, we use separating
channels to balance the accuracy and computational resources. Specifically,
the feature map {${\mathcal{F}_{i}|i\in(2,3,4)}$} with 64 channels is split
evenly into two feature map groups
$\mathcal{F}_{i}^{1},\mathcal{F}_{i}^{2}(i\in(2,3,4))$ and separately fed into
Global Spatial attention (GSA) module and Local Spatial attention (LSA)
module. The outputs of those two attention units are finally concatenated
following by a 1 $\times$ 1 convolution layer. We formulate such a process as

$\displaystyle\mathcal{F}_{i}^{1},\mathcal{F}_{i}^{2}=\rm
Split(\mathcal{F}_{i})$ (4)
$\displaystyle\mathcal{F}_{i}^{{}^{\prime}}=C_{1\times 1}({\rm
Concat}(G_{sa}(\mathcal{F}_{i}^{1}),L_{sa}(\mathcal{F}_{i}^{2}))).$ (5)

Table 1: Quantitative comparison of different methods on Kvasir, ClinicDB, ISIC-2018 and 2018-DSB datasets (seen datasets) to validate our model’s learning ability. $\uparrow$ denotes higher the better and $\downarrow$ denotes lower the better. Methods | Kvasir | ClinicDB | ISIC-2018 | 2018-DSB
---|---|---|---|---
mDice$\uparrow$ | mIou$\uparrow$ | MAE$\downarrow$ | mDice$\uparrow$ | mIou$\uparrow$ | MAE$\downarrow$ | mDice$\uparrow$ | mIou$\uparrow$ | MAE$\downarrow$ | mDice$\uparrow$ | mIou$\uparrow$ | MAE$\downarrow$
U-Net [36] | 0.818 | 0.746 | 0.055 | 0.823 | 0.755 | 0.019 | 0.855 | 0.785 | 0.045 | 0.908 | 0.831 | 0.040
UNet++ [23] | 0.821 | 0.743 | 0.048 | 0.794 | 0.729 | 0.022 | 0.809 | 0.729 | 0.041 | 0.911 | 0.837 | 0.039
PraNet [16] | 0.898 | 0.840 | 0.030 | 0.899 | 0.849 | 0.009 | 0.875 | 0.787 | 0.037 | 0.912 | 0.838 | 0.036
CaraNet [30] | 0.918 | 0.865 | 0.023 | 0.936 | 0.887 | 0.007 | 0.870 | 0.782 | 0.038 | 0.910 | 0.835 | 0.037
TransUNet [10] | 0.913 | 0.857 | 0.028 | 0.935 | 0.887 | 0.008 | 0.880 | 0.809 | 0.036 | 0.915 | 0.845 | 0.033
TransFuse [55] | 0.920 | 0.870 | 0.023 | 0.942 | 0.897 | 0.007 | 0.901 | 0.840 | 0.035 | 0.916 | 0.855 | 0.033
UCTransNet [44] | 0.918 | 0.860 | 0.023 | 0.933 | 0.860 | 0.008 | 0.905 | 0.83 | 0.035 | 0.911 | 0.835 | 0.035
Polyp-PVT [14] | 0.917 | 0.864 | 0.023 | 0.937 | 0.889 | 0.006 | 0.913 | 0.852 | 0.032 | 0.917 | 0.859 | 0.030
DuAT (Ours) | 0.924 | 0.876 | 0.023 | 0.948 | 0.906 | 0.006 | 0.923 | 0.867 | 0.029 | 0.926 | 0.870 | 0.027
Table 2: Quantitative comparison of different methods on ColonDB, ETIS and EndoScene datasets (unseen datasets) to validate our model’s generalization capability. $\uparrow$ denotes higher the better and $\downarrow$ denotes lower the better. Methods | ColonDB | ETIS | EndoScene
---|---|---|---
mDice$\uparrow$ | mIou$\uparrow$ | MAE$\downarrow$ | mDice$\uparrow$ | mIou$\uparrow$ | MAE$\downarrow$ | mDice$\uparrow$ | mIou$\uparrow$ | MAE$\downarrow$
U-Net [36] | 0.512 | 0.444 | 0.061 | 0.398 | 0.335 | 0.036 | 0.710 | 0.627 | 0.022
UNet++ [23] | 0.483 | 0.410 | 0.064 | 0.401 | 0.344 | 0.035 | 0.707 | 0.624 | 0.018
PraNet [16] | 0.712 | 0.640 | 0.043 | 0.628 | 0.567 | 0.031 | 0.851 | 0.797 | 0.010
CaraNet [30] | 0.773 | 0.689 | 0.042 | 0.747 | 0.672 | 0.017 | 0.903 | 0.838 | 0.007
TransUNet [10] | 0.781 | 0.699 | 0.036 | 0.731 | 0.824 | 0.021 | 0.893 | 0.660 | 0.009
TransFuse [55] | 0.781 | 0.706 | 0.035 | 0.737 | 0.826 | 0.020 | 0.894 | 0.654 | 0.009
SSformer [45] | 0.772 | 0.697 | 0.036 | 0.767 | 0.698 | 0.016 | 0.887 | 0.821 | 0.007
Polyp-PVT [14] | 0.808 | 0.727 | 0.031 | 0.787 | 0.706 | 0.013 | 0.900 | 0.833 | 0.007
DuAT (Ours) | 0.819 | 0.737 | 0.026 | 0.822 | 0.746 | 0.013 | 0.901 | 0.840 | 0.005

where $G_{sa}$ denotes the global spatial attention and $L_{sa}$ denotes the
local spatial attention.
$\mathcal{F}_{i}^{{}^{\prime}}\in\mathbb{R}^{\frac{H}{8}\times\frac{W}{8}\times
32}$ is the output features. We will introduce LSA and GSA module in detail in
the following.

$\left(1\right)$ GSA module: The GSA emphasizes the long-range relationship of
each pixel in the spatial space and can be used as a supplement to local
spatial attention. Many efforts [2], [7] claim that the long-range interaction
can make the feature more powerful. Inspired by the manners of extracting
long-range interaction in [2], we simply generate global spatial attention map
($G_{sa}\in\mathbb{R}^{\frac{H}{8}\times\frac{W}{8}\times 32}$) and
$\mathcal{F}_{i}^{1}$ as input as following:

$\displaystyle Att_{G}(\mathcal{F}_{i}^{1})=Softmax(Transpose(C_{1\times
1}(\mathcal{F}_{i}^{1}))),$ (6) $\displaystyle
G_{sa}(\mathcal{F}_{i}^{1})=MLP(Att_{G}(\mathcal{F}_{i}^{1})\otimes\mathcal{F}_{i}^{1})+\mathcal{F}_{i}^{1}.$
(7)

where $Att_{G}(\cdot)$ is the attention operation, $C_{1\times 1}$ means
$1\times 1$ convolution. $\otimes$ denotes matrix multiplication. $MLP(\cdot)$
consists of two fully-connection layers with a ReLU non-linearity and
normalization layer. The first layer of MLP transforms its input to a higher-
dimensional space which the expansion ratio is two, while the second layer
restores the dimension to be the same as the input.

$\left(2\right)$ LSA module: The LSA module extracts the local features of the
region of interest effectively in the spatial dimension of the given feature
map, such as small objects. In short, we compute local spatial attention
response ($L_{sa}\in\mathbb{R}^{\frac{H}{8}\times\frac{W}{8}\times 32}$) and
$\mathcal{F}_{i}^{2}$ as input as follow:

$\displaystyle Att_{L}(\mathcal{F}_{i}^{2})=\sigma(C_{1\times
1}(\mathcal{F}_{c}(\mathcal{F}_{i}^{2}))+\mathcal{F}_{i}^{2})),$ (8)
$\displaystyle
L_{sa}=Att_{L}(\mathcal{F}_{i}^{2})\odot\mathcal{F}_{i}^{2}+\mathcal{F}_{i}^{2}.$
(9)

where $\mathcal{F}_{c}(\cdot)$ denotes cascading three $1\times 1$ convolution
layers and $3\times 3$ depth-wise convolution layers. The number of channels
is adjusted to 32 in the $\mathcal{F}_{c}$. $Att_{L}(\cdot)$ is the local
attention operation, $\sigma(\cdot)$ is the sigmoid function, $\odot$ is
point-wise multiplication. This structural design can efficiently aggregate
local spatial information using fewer parameters.

###  3.4 Loss function

[34, 49] report that combining multiple loss functions with adaptive weights
at different levels can improve the performance of the network with better
convergence speed. Therefore, we use binary cross-entropy loss
($\mathcal{L}^{\omega}_{BCE}(\cdot)$) and the weighted IoU loss (
$\mathcal{L}^{\omega}_{Iou}(\cdot)$) for supervision. Our loss function is
formulated in Eq.10, where $S$ is the two side-outputs ($i,e.,S_{1},S_{2}$)
and $G$ is the ground truth, respectively. $\lambda_{1}$ and $\lambda_{2}$ are
the weighting coefficients.

$\displaystyle\mathcal{L}(S,G)=\lambda_{1}\mathcal{L}^{\omega}_{IoU}(S,G)+\lambda_{2}\mathcal{L}^{\omega}_{BCE}(S,G)$
(10)

Therefore, the total loss $\mathcal{L}_{total}$ for the proposed DuAT can be
formulated as:

$\displaystyle\mathcal{L}_{total}=\mathcal{L}(S_{1},G)+\mathcal{L}(S_{2},G).$
(11)

##  4 Experiments

###  4.1 Datasets

In the experiment, we evaluate our proposed model on three different kinds of
medical image sets: colonoscopy (Colon) images, dermoscopic (Derm) images, and
microscopy (Micro) images, so as to assess the learning ability and
generalization capability of our model. The detailed statistics of each
dataset are shown in Table 3.

Colonoscopy polyp images: Experiments are conducted on five polyp segmentation
datasets (ETIS [43], CVC-ClinicDB (ClinicDB) [37], CVC-ColonDB (ColonDB) [3],
EndoScene-CVC300 (EndoScene) [38], Kvasir-SEG (Kvasir) [22]). We follow the
same training/testing protocols in [14, 16], i.e., the images from the Kvasir
and ClinicDB are randomly split into 80% for training, 10% for validation, and
10% for testing (seen data). And test on the out-of-distribution datasets
which are ColonDB with 380 images, EndoScene with 60 images and ETIS with 196
images (unseen data). Since the resolutions of images are not uniform, we
resize them to 352$\times$352 resolution.

ISIC-2018 Dataset: The dataset comes from ISIC-2018 challenge [13] [41] and is
useful for skin lesion analysis. It includes 2596 images and the corresponding
annotations, which are resized to $512\times 384$ resolution. The images are
randomly split into 80% for training, 10% for validation, and 10% for testing.

2018 Data Science Bowl (2018-DSB): The dataset comes from 2018 Data Science
Bowl challenge [6] and is used to find the nuclei in divergent images,
including 670 images and the corresponding annotations, which are resized to
$256\times 256$ resolution. The images are randomly split into 80% for
training, 10% for validation, and 10% for testing.

Table 3: Statistics on polyp, ISIC-2018 and 2018-DSB datasets. Dataset | Imaging | Images | Shape | Train | Valid | Test
---|---|---|---|---|---|---
Kvasir | Colon | 1000 | Variable | 800 | 100 | 100
ClinicDB | Colon | 612 |  384 $\times$ 288 | 488 | 62 | 62
ColonDB | Colon | 380 |  574 $\times$ 500 | - | - | 380
ETIS | Colon | 196 |  1225 $\times$ 966 | - | - | 196
EndoScene | Colon | 60 |  574 $\times$ 500 | - | - | 60
ISIC-2018 | Derm | 2596 |  512 $\times$ 384 | 2078 | 259 | 259
2018-DSB | Micro | 670 |  256 $\times$ 256 | 536 | 67 | 67

Figure 5: Performance vs. Size on the five polyp datasets. The x-axis is the
proportion size % of polyp and y-axis is the averaged mDice coefficient. Blue
is for our DuAT, orange is for the Polyp-PVT, green is for the TransUnet.

###  4.2 Evaluation Metrics and Implementation Details

Evaluation Metrics. We employ three widely-used metrics i.e., mean Dice
(mDice), mean IoU (mIoU) and mean absolute error (MAE) to evaluate the model
performances. Mean Dice and IoU are the most commonly used metrics and mainly
emphasise the internal consistency of segmentation results. MAE is used to
evaluate the pixel-level accuracy representing the average absolute error
between the prediction and true values.

Implementation Details. We use rotation and horizontal flip for data
augmentation. Considering the differences in the sizes and color of each polyp
image, we adopt a multi-scale training [16, 20] and the color exchange [48].
The network is trained end-to-end by AdamW [29] optimizer. The learning rate
is set to 1e-4 and the weight decay is adjusted to 1e-4 too. The batch size is
set at 16. We use PyTorch framework for implementation with an NVIDIA RTX 3090
GPU. We will provide the source code after the paper is published.

###  4.3 Results

Learning Ability. We first evaluate our proposed DuAT model for its
segmentation performance on seen datasets. As summarized in Table 1, our model
is compared to six recently published models: U-Net [36], UNet++ [23], PraNet
[16], CaraNet [30], TransUNet [10], TransFuse [55], TransFuse [55] and Polyp-
PVT [14] . It can be observed that our DuAT model outperforms all other
models, and achieving 0.924 and 0.948 mean Dice scores on Kvasir and ClinicDB
segmentation respectively. For ISIC-2018 dataset, our DuAT model achieves a
1.0% improvement in terms of mDice and 1.5% of mIoU over SOTA method. For
2018-DSB, DuAT achieves a mIoU of 0.87, mDice of 0.926 and 0.027 of MAE, which
are 1.1%, 1.0%, 0.03% higher than the best performing Polyp-PVT. These results
demonstrate that our model can effectively segment polyps.

Generalization Capabilities. We further evaluate the generalisation capability
of our model on unseen datasets (ETIS, ColonDB, EndoScene). These three
datasets have their own specific challenges and properties. For example,
ColonDB is a small-scale database that contains 380 images from 15 short
colonoscopy sequences. ETIS consists of 196 polyp images for early diagnosis
of colorectal cancer. EndoScene is a re-annotated branch with associated polyp
and background (mucosa and lumen). As seen in Table 2, our model outperforms
the existing medical segmentation baselines on all unseen datasets for all
metrics. Moreover, our DuAT is able to achieve an average dice of 82.2 % on
the most challenging ETIS dataset, 3.5% higher than Polyp-PVT.

Visual Results. We also demonstrate qualitative the performance of our model
on five benchmarks, as given in Figure 4. On ETIS (the first and second row),
DuAT is able to accurately capture the target object’s boundary and detect a
small polyp while other methods fail to detect. On ISIC-2018 (third row), all
methods are able to segment the lesion skin, but our method show the most
similar results compared to the ground truth. On 2018-DSB (the fourth row), we
can observe that our DuAT is better able to capture the presence of nuclei and
obtain better segmentation predictions. More qualitative results can be found
in the supplementary material.

Computational Efficiency. Table LABEL:table122 presents the number of
parameters and floating-point operations for different methods. As our
proposed DuAT and Polyp-PVT [14] adopt the same backbone, they have similar
model size (Params). DuAT uses 24.92M of parameter and 9.88G of FLOPs, which
is more lightweight and compact than CNN-based neural network and Transformer-
based methods.

Table 4: Number of Parameters and FLOPs on a 352 × 352 input image. Note that “N/A” means the result is not available. Method | Type | Params(M) | FLOPs(G)
---|---|---|---
U-Net [36] | CNN | 43.93 | 43.47
UNet++ [23] | CNN | 9.04 | 64.03
PraNet [16] | CNN | 30.50 | 13.01
CaraNet [30] | CNN | 44.59 | 21.65
TransUNet [10] | Transformer | 93.19 | 60.84
SSformer [45] | Transformer | 26.70 | 28.07
TransFuse [55] | Transformer | N/A | N/A
Polyp-PVT [14] | Transformer | 25.08 | 10.00
Ours | Transformer | 24.92 | 9.88
Table 5: Ablation study for DuAT on the Kvasir, ETIS and ISIC-2018 datasets. $\uparrow$ denotes higher the better and $\downarrow$ denotes lower the better. Methods | Kvasir-SEG (seen) | ETIS (unseen) | ISIC-2018 (seen)
---|---|---|---
mDice$\uparrow$ | mIou$\uparrow$ | MAE$\downarrow$ | mDice$\uparrow$ | mIou$\uparrow$ | MAE$\downarrow$ | mDice$\uparrow$ | mIou$\uparrow$ | MAE$\downarrow$
Baseline | 0.910 | 0.856 | 0.030 | 0.759 | 0.668 | 0.035 | 0.877 | 0.783 | 0.040
\+ GSA | 0.912 | 0.860 | 0.029 | 0.772 | 0.675 | 0.030 | 0.887 | 0.803 | 0.038
\+ LSA | 0.916 | 0.863 | 0.028 | 0.785 | 0.690 | 0.027 | 0.900 | 0.839 | 0.035
\+ GSA + LSA (Serial) | 0.914 | 0.863 | 0.028 | 0.786 | 0.695 | 0.025 | 0.909 | 0.845 | 0.034
\+ LSA + GSA (Serial) | 0.910 | 0.860 | 0.029 | 0.799 | 0.713 | 0.021 | 0.910 | 0.852 | 0.033
+GLSA | 0.917 | 0.864 | 0.025 | 0.814 | 0.723 | 0.016 | 0.916 | 0.816 | 0.031
w/o SBA | 0.917 | 0.864 | 0.025 | 0.814 | 0.723 | 0.016 | 0.916 | 0.816 | 0.031
w/o GLSA | 0.915 | 0.863 | 0.026 | 0.790 | 0.696 | 0.023 | 0.901 | 0.800 | 0.033
SBA + GLSA (Ours) | 0.924 | 0.876 | 0.023 | 0.822 | 0.746 | 0.013 | 0.923 | 0.867 | 0.029

Small Object Segmentation Analysis. To demonstrate the detection ability of
our model for small objects, the ratio of the number of pixels in the object
to the number of pixels in the entire image is used to account for the size of
the object. We then evaluate the performance of the segmentation model based
on the size of the object. We set the area with a proportion less than 5%. For
the segmentation model, we first obtain the mean Dice coefficient of the five
polyp datasets. Similar to computing the histogram, we calculate the average
of mean Dice of test data whose size values fall into each interval. For the
small object segmentation analysis, we compare our DuAT with Polyp-PVT and
TransUnet, and the results are given in Figure 5. The overall accuracy of DuAT
is higher than TransUnet [10] and Polyp-PVT [46] on samples with small size
polyps.

###  4.4 Ablation Study

We further conduct ablation study to demonstrate the necessity and
effectiveness of each component of our proposed model on three datasets, and
we choose mDice, mIoU and MAE for evaluation.

Effectiveness of SBA and GLSA. We conduct an experiment to evaluate DuAT
without SBA module “(w/o SBA)”. The performance without the SBA drops sharply
on all three datasets are shown in Table 5. In particular, the mDice is
reduced from 0.822 to 0.814 on ETIS. Moreover, we further investigate the
contribution of the Global-to-Local Spatial Aggregation by removing it from
the overall DuAT and replacing it with convolution operation with a kernel
size of 3, which is denoted as “(w/o GLSA)”. The performance of the complete
DuAT shows an improvement of 2.2 % and 6.7% in terms of mDice and mIoU
respectively on ISIC-2018. After using the two modules (SBA + GLSA), the
model’s performance is improved again. These results demonstrate that these
modules enable our model to distinguish polyp and lesion tissues effectively.

The visual results are given in Figure 6. We observe that the SBA module
facilitates the fine-grained of ambiguous boundaries and the GLSA module
greatly improves the accuracies of small object detection and target object
location.

Figure 6: The effectiveness of each component.

Arrangements of GSA and LSA. GSA and LSA represent global spatial attention
module and local spatial attention module respectively. We further study the
effectiveness and different arrangements of GSA and LSA. The results tested on
Kavsir, ETIS, and ISIC-2018 datasets are shown in Table 5 (the second and
sixth row), and all the methods are using the same backbone PVTv2 [46]. GSA +
LSA(Serial) means first performing GSA then LSA, while LSA + GSA(Serial) is
the opposite. Overall, all improve the baseline and our GLSA group achieves
more accuracy and reliable results. The GLSA module outperforms the GSA, LSA,
GSA + LSA (Serial), LSA + GSA (Serial) by 4.2%, 2.9%, 2.8%, 1.5% in term of
mean Dice on the ETIS dataset.

##  5 Conclusions

In this work, we propose DuAT to address the issues related to medical image
segmentation. Two components, the GLSA and SBA are proposed. Specifically, the
GLSA module extracts the global and local spatial features from the encoder
and is beneficial for locating the large and small objects. The SBA module
alleviates the unclear boundary of high-level features and further improves
its performance. As a result, DuAT can achieve strong learning, generalization
ability, and lightweight segmentation efficiency. Both qualitative and
quantitative results demonstrate the superiority of our DuAT over other
competing methods. We hope that this research will inspire more ideas to solve
the medical image segmentation task and we will extend the proposed model to
tackle 3D medical image segmentation task in the future work.

## Acknowledgments

This work was supported by the Key Program Special Fund in XJTLU (KSF-A-22).

## References

  * [1] Josh Beal, Eric Kim, Eric Tzeng, Dong Huk Park, Andrew Zhai, and Dmitry Kislyuk.  Toward transformer-based object detection.  arXiv preprint arXiv:2012.09958, 2020.
  * [2] Irwan Bello, Barret Zoph, Ashish Vaswani, Jonathon Shlens, and Quoc V Le.  Attention augmented convolutional networks.  In Proceedings of the IEEE/CVF international conference on computer vision, pages 3286–3295, 2019.
  * [3] Jorge Bernal, F Javier Sánchez, Gloria Fernández-Esparrach, Debora Gil, Cristina Rodríguez, and Fernando Vilariño.  Wm-dova maps for accurate polyp highlighting in colonoscopy: Validation vs. saliency maps from physicians.  Computerized Medical Imaging and Graphics, 43:99–111, 2015.
  * [4] Gedas Bertasius, Jianbo Shi, and Lorenzo Torresani.  Semantic segmentation with boundary neural fields.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 3602–3610, 2016.
  * [5] Srinadh Bhojanapalli, Ayan Chakrabarti, Daniel Glasner, Daliang Li, Thomas Unterthiner, and Andreas Veit.  Understanding robustness of transformers for image classification.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 10231–10241, 2021.
  * [6] Juan C Caicedo, Allen Goodman, Kyle W Karhohs, Beth A Cimini, Jeanelle Ackerman, Marzieh Haghighi, CherKeng Heng, Tim Becker, Minh Doan, Claire McQuin, et al.  Nucleus segmentation across imaging experiments: the 2018 data science bowl.  Nature methods, 16(12):1247–1253, 2019.
  * [7] Yue Cao, Jiarui Xu, Stephen Lin, Fangyun Wei, and Han Hu.  Gcnet: Non-local networks meet squeeze-excitation networks and beyond.  In Proceedings of the IEEE/CVF International Conference on Computer Vision Workshops, pages 0–0, 2019.
  * [8] Nicolas Carion, Francisco Massa, Gabriel Synnaeve, Nicolas Usunier, Alexander Kirillov, and Sergey Zagoruyko.  End-to-end object detection with transformers.  In European conference on computer vision, pages 213–229. Springer, 2020.
  * [9] Hanting Chen, Yunhe Wang, Tianyu Guo, Chang Xu, Yiping Deng, Zhenhua Liu, Siwei Ma, Chunjing Xu, Chao Xu, and Wen Gao.  Pre-trained image processing transformer.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 12299–12310, 2021.
  * [10] Jieneng Chen, Yongyi Lu, Qihang Yu, Xiangde Luo, Ehsan Adeli, Yan Wang, Le Lu, Alan L Yuille, and Yuyin Zhou.  Transunet: Transformers make strong encoders for medical image segmentation.  arXiv preprint arXiv:2102.04306, 2021.
  * [11] Liang-Chieh Chen, George Papandreou, Iasonas Kokkinos, Kevin Murphy, and Alan L Yuille.  Deeplab: Semantic image segmentation with deep convolutional nets, atrous convolution, and fully connected crfs.  IEEE transactions on pattern analysis and machine intelligence, 40(4):834–848, 2017.
  * [12] Liang-Chieh Chen, Yukun Zhu, George Papandreou, Florian Schroff, and Hartwig Adam.  Encoder-decoder with atrous separable convolution for semantic image segmentation.  In Proceedings of the European conference on computer vision (ECCV), pages 801–818, 2018.
  * [13] Noel Codella, Veronica Rotemberg, Philipp Tschandl, M Emre Celebi, Stephen Dusza, David Gutman, Brian Helba, Aadi Kalloo, Konstantinos Liopyris, Michael Marchetti, et al.  Skin lesion analysis toward melanoma detection 2018: A challenge hosted by the international skin imaging collaboration (isic).  arXiv preprint arXiv:1902.03368, 2019.
  * [14] Bo Dong, Wenhai Wang, Deng-Ping Fan, Jinpeng Li, Huazhu Fu, and Ling Shao.  Polyp-pvt: Polyp segmentation with pyramid vision transformers.  arXiv preprint arXiv:2108.06932, 2021.
  * [15] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  arXiv preprint arXiv:2010.11929, 2020.
  * [16] Deng-Ping Fan, Ge-Peng Ji, Tao Zhou, Geng Chen, Huazhu Fu, Jianbing Shen, and Ling Shao.  Pranet: Parallel reverse attention network for polyp segmentation.  In International conference on medical image computing and computer-assisted intervention, pages 263–273. Springer, 2020.
  * [17] Pasqualino Favoriti, Gabriele Carbone, Marco Greco, Felice Pirozzi, Raffaele Emmanuele Maria Pirozzi, and Francesco Corcione.  Worldwide burden of colorectal cancer: a review.  Updates in surgery, 68(1):7–11, 2016.
  * [18] Huazhu Fu, Jun Cheng, Yanwu Xu, Damon Wing Kee Wong, Jiang Liu, and Xiaochun Cao.  Joint optic disc and cup segmentation based on multi-label deep network and polar transformation.  IEEE transactions on medical imaging, 37(7):1597–1605, 2018.
  * [19] Kai Han, An Xiao, Enhua Wu, Jianyuan Guo, Chunjing Xu, and Yunhe Wang.  Transformer in transformer.  Advances in Neural Information Processing Systems, 34, 2021.
  * [20] Chien-Hsiang Huang, Hung-Yu Wu, and Youn-Long Lin.  Hardnet-mseg: a simple encoder-decoder polyp segmentation neural network that achieves over 0.9 mean dice and 86 fps.  arXiv preprint arXiv:2101.07172, 2021.
  * [21] Zilong Huang, Xinggang Wang, Lichao Huang, Chang Huang, Yunchao Wei, and Wenyu Liu.  Ccnet: Criss-cross attention for semantic segmentation.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 603–612, 2019.
  * [22] Debesh Jha, Pia H Smedsrud, Michael A Riegler, Pål Halvorsen, Thomas de Lange, Dag Johansen, and Håvard D Johansen.  Kvasir-seg: A segmented polyp dataset.  In International Conference on Multimedia Modeling, pages 451–462. Springer, 2020.
  * [23] Debesh Jha, Pia H Smedsrud, Michael A Riegler, Dag Johansen, Thomas De Lange, Pål Halvorsen, and Håvard D Johansen.  Resunet++: An advanced architecture for medical image segmentation.  In 2019 IEEE International Symposium on Multimedia (ISM), pages 225–2255. IEEE, 2019.
  * [24] Ge-Peng Ji, Lei Zhu, Mingchen Zhuge, and Keren Fu.  Fast camouflaged object detection via edge-based reversible re-calibration network.  Pattern Recognition, 123:108414, 2022.
  * [25] Xiangtai Li, Xia Li, Li Zhang, Guangliang Cheng, Jianping Shi, Zhouchen Lin, Shaohua Tan, and Yunhai Tong.  Improving semantic segmentation via decoupled body and edge supervision.  In European Conference on Computer Vision, pages 435–452. Springer, 2020.
  * [26] Xiangtai Li, Li Zhang, Ansheng You, Maoke Yang, Kuiyuan Yang, and Yunhai Tong.  Global aggregation then local distribution in fully convolutional networks.  arXiv preprint arXiv:1909.07229, 2019.
  * [27] Xiangtai Li, Houlong Zhao, Lei Han, Yunhai Tong, Shaohua Tan, and Kuiyuan Yang.  Gated fully fusion for semantic segmentation.  In Proceedings of the AAAI conference on artificial intelligence, volume 34, pages 11418–11425, 2020.
  * [28] Ze Liu, Yutong Lin, Yue Cao, Han Hu, Yixuan Wei, Zheng Zhang, Stephen Lin, and Baining Guo.  Swin transformer: Hierarchical vision transformer using shifted windows.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 10012–10022, 2021.
  * [29] Ilya Loshchilov and Frank Hutter.  Decoupled weight decay regularization.  arXiv preprint arXiv:1711.05101, 2017.
  * [30] Ange Lou, Shuyue Guan, and Murray Loew.  Caranet: Context axial reverse attention network for segmentation of small medical objects.  arXiv preprint arXiv:2108.07368, 2021.
  * [31] Yi Lu, Yaran Chen, Dongbin Zhao, and Jianxin Chen.  Graph-fcn for image semantic segmentation.  In International symposium on neural networks, pages 97–105. Springer, 2019.
  * [32] Haoxiang Ma, Hongyu Yang, and Di Huang.  Boundary guided context aggregation for semantic segmentation.  arXiv preprint arXiv:2110.14587, 2021.
  * [33] Prashant Mathur, Krishnan Sathishkumar, Meesha Chaturvedi, Priyanka Das, Kondalli Lakshminarayana Sudarshan, Stephen Santhappan, Vinodh Nallasamy, Anish John, Sandeep Narasimhan, Francis Selvaraj Roselind, et al.  Cancer statistics, 2020: report from national cancer registry programme, india.  JCO Global Oncology, 6:1063–1075, 2020.
  * [34] Xuebin Qin, Zichen Zhang, Chenyang Huang, Chao Gao, Masood Dehghan, and Martin Jagersand.  Basnet: Boundary-aware salient object detection.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 7479–7489, 2019.
  * [35] René Ranftl, Alexey Bochkovskiy, and Vladlen Koltun.  Vision transformers for dense prediction.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 12179–12188, 2021.
  * [36] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In International Conference on Medical image computing and computer-assisted intervention, pages 234–241. Springer, 2015.
  * [37] Juan Silva, Aymeric Histace, Olivier Romain, Xavier Dray, and Bertrand Granado.  Toward embedded detection of polyps in wce images for early diagnosis of colorectal cancer.  International journal of computer assisted radiology and surgery, 9(2):283–293, 2014.
  * [38] Nima Tajbakhsh, Suryakanth R Gurudu, and Jianming Liang.  Automated polyp detection in colonoscopy videos using shape and context information.  IEEE transactions on medical imaging, 35(2):630–644, 2015.
  * [39] Towaki Takikawa, David Acuna, Varun Jampani, and Sanja Fidler.  Gated-scnn: Gated shape cnns for semantic segmentation.  In Proceedings of the IEEE/CVF international conference on computer vision, pages 5229–5238, 2019.
  * [40] Hugo Touvron, Matthieu Cord, Matthijs Douze, Francisco Massa, Alexandre Sablayrolles, and Hervé Jégou.  Training data-efficient image transformers & distillation through attention.  In International Conference on Machine Learning, pages 10347–10357. PMLR, 2021.
  * [41] Philipp Tschandl, Cliff Rosendahl, and Harald Kittler.  The ham10000 dataset, a large collection of multi-source dermatoscopic images of common pigmented skin lesions.  Scientific data, 5(1):1–9, 2018.
  * [42] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  Advances in neural information processing systems, 30, 2017.
  * [43] David Vázquez, Jorge Bernal, F Javier Sánchez, Gloria Fernández-Esparrach, Antonio M López, Adriana Romero, Michal Drozdzal, and Aaron Courville.  A benchmark for endoluminal scene segmentation of colonoscopy images.  Journal of healthcare engineering, 2017, 2017.
  * [44] Haonan Wang, Peng Cao, Jiaqi Wang, and Osmar R Zaiane.  Uctransnet: rethinking the skip connections in u-net from a channel-wise perspective with transformer.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pages 2441–2449, 2022.
  * [45] Jinfeng Wang, Qiming Huang, Feilong Tang, Jia Meng, Jionglong Su, and Sifan Song.  Stepwise feature fusion: Local guides global.  arXiv preprint arXiv:2203.03635, 2022.
  * [46] Wenhai Wang, Enze Xie, Xiang Li, Deng-Ping Fan, Kaitao Song, Ding Liang, Tong Lu, Ping Luo, and Ling Shao.  Pvt v2: Improved baselines with pyramid vision transformer.  Computational Visual Media, pages 1–10, 2022.
  * [47] Xiaolong Wang, Ross Girshick, Abhinav Gupta, and Kaiming He.  Non-local neural networks.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 7794–7803, 2018.
  * [48] Jun Wei, Yiwen Hu, Ruimao Zhang, Zhen Li, S Kevin Zhou, and Shuguang Cui.  Shallow attention network for polyp segmentation.  In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 699–708. Springer, 2021.
  * [49] Jun Wei, Shuhui Wang, and Qingming Huang.  F3net: fusion, feedback and focus for salient object detection.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 12321–12328, 2020.
  * [50] Enze Xie, Wenhai Wang, Zhiding Yu, Anima Anandkumar, Jose M Alvarez, and Ping Luo.  Segformer: Simple and efficient design for semantic segmentation with transformers.  Advances in Neural Information Processing Systems, 34, 2021.
  * [51] Fisher Yu and Vladlen Koltun.  Multi-scale context aggregation by dilated convolutions.  arXiv preprint arXiv:1511.07122, 2015.
  * [52] Li Zhang, Xiangtai Li, Anurag Arnab, Kuiyuan Yang, Yunhai Tong, and Philip HS Torr.  Dual graph convolutional network for semantic segmentation.  arXiv preprint arXiv:1909.06121, 2019.
  * [53] Li Zhang, Dan Xu, Anurag Arnab, and Philip HS Torr.  Dynamic graph message passing networks.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 3726–3735, 2020.
  * [54] Ruifei Zhang, Guanbin Li, Zhen Li, Shuguang Cui, Dahong Qian, and Yizhou Yu.  Adaptive context selection for polyp segmentation.  In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 253–262. Springer, 2020.
  * [55] Yundong Zhang, Huiye Liu, and Qiang Hu.  Transfuse: Fusing transformers and cnns for medical image segmentation.  In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 14–24. Springer, 2021.
  * [56] Zhenli Zhang, Xiangyu Zhang, Chao Peng, Xiangyang Xue, and Jian Sun.  Exfuse: Enhancing feature fusion for semantic segmentation.  In Proceedings of the European conference on computer vision (ECCV), pages 269–284, 2018.
  * [57] Mingmin Zhen, Jinglu Wang, Lei Zhou, Shiwei Li, Tianwei Shen, Jiaxiang Shang, Tian Fang, and Long Quan.  Joint semantic segmentation and boundary detection using iterative pyramid contexts.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 13666–13675, 2020.
  * [58] Sixiao Zheng, Jiachen Lu, Hengshuang Zhao, Xiatian Zhu, Zekun Luo, Yabiao Wang, Yanwei Fu, Jianfeng Feng, Tao Xiang, Philip HS Torr, et al.  Rethinking semantic segmentation from a sequence-to-sequence perspective with transformers.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 6881–6890, 2021.
  * [59] Zongwei Zhou, Md Mahfuzur Rahman Siddiquee, Nima Tajbakhsh, and Jianming Liang.  Unet++: A nested u-net architecture for medical image segmentation.  In Deep learning in medical image analysis and multimodal learning for clinical decision support, pages 3–11. Springer, 2018.


# Numerical and modeling error assessment of large-eddy simulation using
direct-numerical-simulation-aided large-eddy simulation

H. Jane Bae<EMAIL_ADDRESS>Adrián Lozano-Durán<EMAIL_ADDRESS>Graduate Aerospace Laboratories, California Institute of Technology, Pasadena,
California 91125, USA  Department of Aeronautics and Astronautics,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA

###### Abstract

We study the numerical errors of large-eddy simulation (LES) in isotropic and
wall-bounded turbulence. A direct-numerical-simulation (DNS)-aided LES
formulation, where the subgrid-scale (SGS) term of the LES is computed by
using filtered DNS data is introduced. We first verify that this formulation
has zero error in the absence of commutation error between the filter and the
differentiation operator of the numerical algorithm. This method allows the
evaluation of the time evolution of numerical errors for various numerical
schemes at grid resolutions relevant to LES. The analysis shows that the
numerical errors are of the same order of magnitude as the modeling errors and
often cancel each other. This supports the idea that supervised machine
learning algorithms trained on filtered DNS data might not be suitable for
robust SGS model development, as this approach disregards the existence of
numerical errors in the system that accumulates over time. The assessment of
errors in turbulent channel flow also identifies that numerical errors close
to the wall dominate, which has implications for the development of wall
models.

##  1 Introduction

DNS of most turbulent flows relevant for industrial applications is not
tractable because the range of scales of motions in these flows is so large
that the computational cost becomes prohibitive. In LES, the effect of the
small scales on the larger ones is modeled through an SGS model, reducing the
computational cost by several orders of magnitude. The accuracy of the
representation of turbulent flows depends on both the performance or modeling
capabilities of the choice of the SGS model, the adequacy of the LES
resolution, and the numerical method employed to discretize the equations. In
addition, the various sources of error can interact, which may lead to further
complications in the validation of SGS models in LES. Given that accurate
predictions are required in many engineering and scientific applications, a
close assessment of LES numerical errors is crucial. This will also have
implications for developing SGS models using supervised learning methods,
which most often use filtered DNS data to train LES models without regarding
the numerical errors.

Various SGS models have been introduced in the past decades such as the eddy-
viscosity models [1, 2, 3, 4, 5], stochastic models [6, 7], similarity models
[8, 9, 10], and mixed models [11, 12], among others. In addition, an important
development in SGS modeling arose with the dynamic procedure [13, 14, 15],
which allows optimization of parameters in the eddy-viscosity SGS models in
accordance with the turbulent flow that is simulated. First works aiming to
assess the accuracy of SGS models include the pioneering investigation by
Clark et al. [8], who established the numerical study of decaying isotropic
turbulence as a reference benchmark, although the grid resolutions and
Reynolds numbers tested were highly constrained by the computational resources
of the time. Since then, common benchmarks for LES have broadened to include
simple hydrodynamic cases such as forced or decaying isotropic turbulence [4],
rotating homogeneous turbulence [16], spatial or temporal mixing layers [17,
18] and plane turbulent channel flow [19, 13, 20, 21], among others. See
Bonnet et al. [22] for an overview of cases for LES validation.

The validation cases mentioned above are representative of most LES error
quantification, where the numerical error and the modeling error are not
distinguished. Vreman et al. [23] was one of the first to quantify the
modeling and discretization errors separately for various flow properties. In
this analysis, modeling and numerical errors were found to be of comparable
magnitude and could partially cancel each other. Further error analysis for
LES was presented in Geurts and Fröhlich [24] in the context of the “subgrid-
activity” parameter. More recently, Meyers et al. [25] differentiated the two
errors and again showed that the partial cancellation of both sources can lead
to coincidental accurate results. Along the same line, Meyers and Sagaut [26]
studied the combined effect of discretization and model errors, and a further
series of works resulted in the error-landscape-methodology framework reviewed
by Meyers [27], where it is stressed that the determination of LES quality
based on one single metric alone may provide misleading results. However, in
these works, the LES errors are quantified as the error occurring in SGS
models on a DNS grid by setting the filter size of the SGS grid to the LES
grid size. As most SGS models for LES are not designed or expected to work in
DNS grid resolutions, the error defined above may not be an accurate depiction
of the true modeling error observed in LES grid resolutions. Moreover, this
method is able to assess SGS models when there is a clearly defined filter
size, which is not straightforward for anisotropic and unstructured grids
[28]. Ideally, the modeling error would compare a filtered DNS solution in an
LES grid with the LES solution. This can be made possible by first quantifying
the correct numerical error for an LES using explicitly filtered LES.

Previous works on explicitly filtered LES include the study of Winckelmans et
al. [29], who investigated a two-dimensional explicitly filtered isotropic
turbulence and channel flow LES to evaluate various mixed subgrid/subfilter
scale models. Stolz et al. [30] implemented the three-dimensional filtering
schemes of Vasilyev et al. [31] by using an approximate deconvolution model
for the convective terms in the LES equations. Lund [32] applied two-
dimensional explicit filters to a channel flow and evaluated the performance
of explicitly filtered versus implicitly filtered LES. Gullbrand and Chow [33]
attempted the first grid-independent solution of the LES equations with
explicit filtering. Bose et al. [34] further investigated the grid
independence of explicitly filtered LES with a three-dimensional filter for
turbulent channel flows. The recent works of Bae and Lozano-Durán [35, 36]
introduce a DNS-aided LES (DAL) framework using explicitly filtered LES that
computes the exact SGS stress term in the absence of commutation errors from
filtered DNS data.

From a different perspective, machine-learning-based SGS models trained using
high-fidelity data, typically DNS, have recently emerged (see reviews [37, 38,
39, 40]). SGS models have been developed in several canonical and complex
flows, including two-dimensional homogeneous isotropic turbulence [41, 42],
three-dimensional forced and [43, 44, 45, 46] decaying [47, 48] isotropic
turbulence, turbulent channel flow [49, 50, 51], and flow over an aircraft
[52]. A majority of the research utilizes a supervised learning framework,
where the model is trained to produce accurate SGS stresses obtained from DNS
(or filtered DNS) data, and is often based on single-step target values to
limit computational challenges in minimizing the model prediction error.
Therefore, it is important to differentiate between _a priori_ and _a
posteriori_ testing for these models. While the model might perform well in _a
priori_ testing, _a posteriori_ testing is performed by integrating the LES
equations in time. Due to the single-step cost function as well as the
discrepancy in the filtered Navier-Stokes and LES equations, the resulting
machine-learned model is not trained to compensate for the discrepancies
between DNS and LES and the compounding (numerical and modeling) errors. The
issue of inconsistency in data-driven SGS models has been exposed by studies
that performed _a posteriori_ testing [53, 49, 48]. Thus, in order to account
for the discrepancies between DNS and LES in machine-learning-based model
development, a systematic understanding of the different errors associated
with LES SGS modeling is necessary.

It is the aim of this paper to introduce a framework that will allow the
evaluation of the time evolution of numerical errors in LES. For that purpose,
we introduce a method that produces zero modeling error using the DAL
framework [35, 36] as a method to define separate numerical and modeling
errors. Our goal is to assess the errors for unbounded and wall-bounded
turbulent flows using this methodology. In particular, for wall-bounded flows,
the wall-normal dependency of errors in LES will be investigated.

The paper is organized as follows. In Section 2, the filtered Navier-Stokes
equations are compared to the traditional LES equations to clearly define
numerical error, and the DAL framework is introduced. We then introduce the
numerical experiments in Section 3, where we evaluate the numerical errors for
forced isotropic turbulence and turbulent channel flows. Finally, we conclude
the paper in Section 4.

##  2 Zero-modeling-error formulation

###  2.1 Filtered Navier-Stokes equations and large-eddy simulation equations

The incompressible Navier-Stokes equations are given by

${\frac{\partial{u}_{i}}{\partial t}}+{\frac{\partial{u}_{i}{u}_{j}}{\partial
x_{j}}}=-\frac{1}{\rho}{\frac{\partial{p}}{\partial
x_{i}}}+\nu{\frac{\partial^{2}{u}_{i}}{\partial x_{j}\partial
x_{j}}},\quad{\frac{\partial{u}_{i}}{\partial x_{i}}}=0,$ (1)

where the velocity components are represented by ${u}_{i}$ for $i=1,2,3$ (or
equivalently $u$, $v$, and $w$) for the directions $x_{i}$ ($x$, $y$, and
$z$), $\rho$ is the fluid density, $\nu$ is the kinematic viscosity, and $p$
is the pressure. As LES is formally derived from the filtered Navier-Stokes
equations, we define a filter operator on a variable $\phi$ in integral form

$\bar{\phi}(\boldsymbol{x})=\int_{\Omega}G(t,\boldsymbol{x},\boldsymbol{x^{\prime}})\phi(\boldsymbol{x^{\prime}})\mathrm{d}\boldsymbol{x^{\prime}},$
(2)

where $\boldsymbol{x}=(x_{1},x_{2},x_{3})$, $G$ is the filter kernel, and
$\Omega$ is the domain of integration. When Eq. (1) is filtered with Eq. (2),
the resulting equations are

$\overline{\frac{\partial{u}_{i}}{\partial
t}}+\overline{\frac{\partial{u}_{i}{u}_{j}}{\partial
x_{j}}}=-\frac{1}{\rho}\overline{\frac{\partial{p}}{\partial
x_{i}}}+\nu\overline{\frac{\partial^{2}{u}_{i}}{\partial x_{j}\partial
x_{j}}},\quad\overline{\frac{\partial{u}_{i}}{\partial x_{i}}}=0.$ (3)

However, in LES, we have access to the filtered velocity quantities
$\bar{u}_{i}$ and not the full flow field $u_{i}$. In order to be able to
solve Eq. (3) for a practical LES applications, the order of the filter and
the differentiation operator must be inverted such that the equations become

${\frac{\partial\bar{u}_{i}}{\partial
t}}+{\frac{\partial\overline{\bar{u}_{i}\bar{u}_{j}}}{\partial
x_{j}}}=-\frac{1}{\rho}{\frac{\partial\bar{p}}{\partial
x_{i}}}+\nu{\frac{\partial^{2}\bar{u}_{i}}{\partial x_{j}\partial
x_{j}}}-\frac{\partial\mathcal{T}_{ij}}{\partial
x_{j}}+\delta^{\text{num}}_{i},\quad{\frac{\partial\bar{u}_{i}}{\partial
x_{i}}}+\delta^{\text{num}}=0,$ (4)

where $\delta^{\text{num}}_{i}$ and $\delta^{\text{num}}$ are the commutation
error due to the inversion of the order of operators. In traditional LES, the
commutation error is often neglected or assumed to be zero regardless of the
choice of filter and differentiation operator. However, in most numerical
methods used for LES, these commutation errors are not negligible even with
the use of commutative filters [54] as the commutation of the two operators is
guaranteed only up to a finite order. Here, we retain
$\delta^{\text{num}}_{i}$ and $\delta^{\text{num}}$ to signify the difference
between the formulations in Eq. (3) and Eq. (4). The term
$\mathcal{T}_{ij}=\overline{u_{i}u}_{j}-\overline{\bar{u}_{i}\bar{u}_{j}}$ is
the SGS stress tensor, which is modeled through an SGS model in LES. At each
time step, the effect of using an SGS model can be characterized as
$\delta^{\text{mod}}_{i}=\partial/\partial
x_{j}\left(\mathcal{T}_{ij}^{\text{SGS}}-\mathcal{T}_{ij}\right)$, where
$\mathcal{T}_{ij}^{\text{SGS}}$ is the modeled SGS stress tensor given by the
choice of SGS model. In this paper, we will categorize the time-integrated
effect of $\delta^{\text{num}}_{i}$ and $\delta^{\text{num}}$ as the numerical
error, $\mathcal{E}^{\text{num}}$ in the absence of $\delta_{i}^{\text{mod}}$.

###  2.2 Zero-modeling error through DNS-aided LES

Figure 1: Schematic of DAL. The SGS stress tensor
$\mathcal{T}_{ij}^{\text{DAL}}$ is computed at each time step from the DNS
flow field at the same time step and then applied to the LES equations.

One possible way of enforcing zero modeling error in LES is by using an exact
SGS model that can be produced by running a DNS concurrently. That is,
starting from an initial DNS velocity field $u^{\text{DNS}}_{i}(t_{0})$ and
filtering it, we can compute the initial condition
$\bar{u}^{\text{DAL}}_{i}(t_{0})$ for the DAL simulations such that
$\overline{u^{\text{DNS}}_{i}}(t_{0})=\bar{u}^{\text{DAL}}_{i}(t_{0})$. At any
given time step $t_{n}$ of the DAL, the exact SGS stress tensor
$\mathcal{T}_{ij}$ can be given by the DNS at the same time step $t_{n}$,
i.e.,

$\mathcal{T}_{ij}^{\text{DAL}}(t_{n})=\overline{u^{\text{DNS}}_{i}(t_{n})u^{\text{DNS}}_{j}(t_{n})}-\overline{\bar{u}^{\text{DAL}}_{i}(t_{n}){\bar{u}^{\text{DAL}}}_{j}(t_{n})}.$
(5)

A schematic of the proposed method is given in Figure 1. Note that this would
be not possible for traditional implicitly filtered LES due to the lack of an
explicit filter operator. In the absence of numerical (commutation) error, the
filtered DNS and DAL solutions will be identical, thus the name zero modeling
error. However, in the presence of numerical error, the solutions of
$\overline{u^{\text{DNS}}_{i}}$ and $\bar{u}_{i}^{\text{DAL}}$ will diverge
from the integrated effect of $\delta_{i}^{\text{num}}$, and the difference in
the filtered DNS and DAL solution will provide a measure for the numerical
error of the simulation without the presence of modeling error in a grid
relevant to LES.

##  3 Error evaluation in LES

We consider two test cases, forced isotropic turbulence and turbulent channel
flow, to evaluate the numerical and modeling errors of LES. We quantify the
error in the LES flow fields as the $L_{2}$ error of the filtered DNS and LES
velocity fields,

$\displaystyle\mathcal{E}^{\text{tot}}(t)=$
$\displaystyle\bigg{[}\frac{1}{L_{x}L_{y}L_{z}}\iiint\sum_{i=1}^{3}\left(\overline{u_{i}^{\text{DNS}}}(t)-\bar{u}_{i}^{\text{LES}}(t)\right)^{2}\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z\bigg{]}^{1/2},$
(6)

where $L_{x,y,z}$ is the length of the domain in the corresponding direction.
The error $\mathcal{E}^{\text{tot}}$ can be seen as a combination of the
numerical and modeling error by considering the square-root of the integrand
of Eq. (6) as

$\overline{u_{i}^{\text{DNS}}}(t)-\bar{u}_{i}^{\text{LES}}(t)=\left(\overline{u_{i}^{\text{DNS}}}(t)-\bar{u}_{i}^{\text{DAL}}(t)\right)+\Big{(}\bar{u}_{i}^{\text{DAL}}(t)-\bar{u}_{i}^{\text{LES}}(t)\Big{)},$
(7)

where we define the numerical error as

$\mathcal{E}^{\text{num}}(t)=\left[\frac{1}{L_{x}L_{y}L_{z}}\iiint\sum_{i=1}^{3}\left(\overline{u_{i}^{\text{DNS}}}(t)-\bar{u}_{i}^{\text{DAL}}(t)\right)^{2}\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z\right]^{1/2}.$
(8)

In essence, the numerical error quantifies the error due to the filter
operator.

In the case of channel flow, we also quantify a wall-normal distance dependent
error

$\mathcal{E}^{\text{tot}}_{y}(y,t)=\left[\frac{1}{L_{x}L_{z}}\iint\sum_{i=1}^{3}\left(\overline{u_{i}^{\text{DNS}}}(y,t)-\bar{u}_{i}^{\text{LES}}(y,t)\right)^{2}\mathrm{d}x\,\mathrm{d}z\right]^{1/2},$
(9)

where $y$ is the wall-normal direction. An analogous $y$-dependent numerical
error can be defined as $\mathcal{E}^{\text{num}}_{y}$.

It is worth remarking the fact that if an SGS model is able to counteract the
numerical errors with the modeling errors such that
$\left(\overline{u_{i}^{\text{DNS}}}(t)-\bar{u}_{i}^{\text{DAL}}(t)\right)=-\Big{(}\bar{u}_{i}^{\text{DAL}}(t)-\bar{u}_{i}^{\text{LES}}(t)\Big{)}$,
the total error is zero. Thus, not only does the DAL formulation provide a
method to evaluate errors, it can be used to inform future SGS models.
However, SGS models are not expected to provide the correct closure term in an
instantaneous manner. Rather, the goal of LES is to obtain the correct
statistics of the system. Thus, computing instantaneous modeling errors does
not align with the goal of LES and SGS modeling and will not be considered in
this paper.

###  3.1 Forced isotropic turbulence

(a)  (b)

Figure 2: (a) Time evolution of the numerical error $\mathcal{E}^{\text{num}}$
using dealiased Fourier discretization in space and the Fourier cutoff filter
for grid resolutions $32^{3}$ ( ), $64^{3}$ ( ), and $128^{3}$ ( ). Vertical
lines ( ) divide the initial response region and growth region. (b) Time
evolution of the numerical error $\mathcal{E}^{\text{num}}$ using a second-
order finite difference in space and downsampling for grid resolutions
$80^{3}$ ( ) and $48^{3}$ times ( ). Vertical lines ( ) divide the initial
response region, growth region, and saturation region.

To demonstrate the method given in section 2.2 provides zero modeling error
and evaluate modeling error, we first perform a simulation of forced isotropic
turbulence using a dealiased Fourier discretization in space in a triply
periodic domain of size $2\pi^{3}$ and advanced in time using a fourth-order
Runge-Kutta time stepping method. A linear forcing [55], $f_{i}=Au_{i}$, is
applied with $A=0.45$. The DNS is performed using a $256^{3}$ mesh at
$Re_{\lambda}\approx 100$. Three DAL (section 2.2) cases with grid resolution
of $32^{2}$, $64^{3}$ and $128^{3}$ are performed with the filter operator
given by a sharp Fourier cutoff filter at the largest wavenumber resolvable by
the LES grid. The LES simulations are run with the same time-step size as DNS
to take advantage of the DAL formulation without interpolation in time. The
initial condition for the numerical experiment was taken from a separate DNS
after transients. The errors are normalized by the square-root of the
turbulence kinetic energy $k$, and the time is given in integral times
$k/|\varepsilon|$, where $\varepsilon$ is the dissipation rate.

We evaluate the numerical error using DAL. Since DAL can be thought of as
supplying the exact SGS model that produces zero modeling error, the
difference between the filtered DNS flow field and the DAL flow field can be
considered the numerical error as defined in the beginning of section 2.1. The
sharp Fourier cutoff filter used in the simulations commutes with the
derivative operator in Fourier space, which gives, in theory, zero numerical
error. In practice, the numerical error is zero only up to machine precision,
and the resulting small perturbations in the velocity profile will grow
exponentially following the leading Lyapunov exponent of the Navier-Stokes
equation until it reaches saturation [56]. The Lyapunov exponent $\lambda$
represents the rate of separation, such that the error grows proportional to
$\exp(\lambda t)$, and its reciprocal is closely related to the predictability
horizon of a chaotic solution, such as solutions to the Navier-Stokes
equations. Thus, DAL solution should give $\mathcal{E}^{\text{num}}\approx 0$
initially with a small exponential growth rate to ensure that this formulation
is effective.

In Figure 2(a), we plot $\mathcal{E}^{\text{num}}(t)$ for the cases
considered. As expected, the error is initially $O(10^{-15})$ for all LES grid
resolutions, with the error following a linear curve, known as the initial
response, until $t|\varepsilon|/k\approx 1.5$. Then, the error grows
exponentially in time with an exponential growth rate, or the Lyapunov
exponent, of $\lambda k/|\varepsilon|\approx 0.4$. This demonstration shows
that for short times, the two flow fields can be considered identical (the
error will be less than $10^{-10}$ for 30 integral times), and thus, the
modeling error evaluated from DAL flow fields is zero in the absence of
numerical error.

We repeat the previous experiment using a second-order finite difference
staggered discretization in space and advanced in time using a third-order
Runge-Kutta time-stepping method. The same linear forcing is applied using a
$240^{3}$ mesh. The DAL cases are three and five times coarser than the DNS
case in all three spatial directions, i.e., computed on a $80^{3}$ and
$48^{3}$ mesh. Due to the numerical method of the simulation, there is no
filter operator that maps a velocity field of a DNS grid to a smaller LES grid
and commutes with the differentiation matrix (see Appendix A). An example of a
case with the LES grid equal to the DNS grid is given in Bae and Lozano-Durán
[36] for a differential filter [57] utilizing the extension method [35, 36].
However, for practical applications, LES should reduce the computational cost
by using a low pass filter and a smaller computational domain. Here, we
utilize a down-sampling method where the velocity field at a certain point of
the LES is the same as the (interpolated) DNS velocity field at the same
point.

The numerical error $\mathcal{E}^{\text{num}}$ as a function of time is given
in figure 2(b). Unlike the case using Fourier discretization where the
differentiation operator commutes with the Fourier-cutoff filter, the down-
sampling does not commute with the second-order finite-difference operation.
The error is much larger from the first time step, assuming a polynomial
growth with a growth factor of $0.4$ until it reaches saturation. Unlike the
previous case, the initial error is too large and the gap between the initial
and fully saturated error is too small to develop an exponential growth
region. Instead, the polynomial growth rate shows how the initial error
develops.

###  3.2 Channel flow

We perform a set of plane turbulent channel simulations at friction Reynolds
number $Re_{\tau}=u_{\tau}h/\nu\approx 180$ and $550$, where $h$ is the
channel half-height and $u_{\tau}$ is the friction velocity at the wall. The
simulations are computed with a staggered second-order finite difference [58]
and a fractional-step method [59] with a third-order Runge-Kutta time-
advancing scheme [60]. The code has been validated in previous studies in
turbulent channel flows [61, 62, 63]. The size of the channel is $2\pi h\times
2h\times\pi h$ in the $x$, $y$ and $z$ directions, which are the streamwise,
wall-normal, and spanwise directions, respectively. It has been shown that
this domain size is sufficient to accurately predict one-point statistics for
$Re_{\tau}$ up to 4200 [64]. Periodic boundary conditions are imposed in the
streamwise and spanwise directions, and the flow is driven by imposing a
constant mean pressure gradient.

The DNS grid resolution utilizes $192\times 120\times 96$ grid points for the
$Re_{\tau}\approx 180$ case and $510\times 300\times 510$ grid points for the
$Re_{\tau}\approx 550$ case in the three spatial directions. The streamwise
and spanwise directions are uniform with $\Delta_{x}^{+}\approx 6$ and
$\Delta_{z}^{+}\approx 6$ or $3$, where the superscript $+$ denotes wall units
given by $\nu$ and $u_{\tau}$. Non-uniform meshes are used in the normal
direction, with the grid stretched toward the wall according to a hyperbolic
tangent distribution. The height of the first grid cell at the wall is
$\Delta_{y}^{+}\approx 0.2$ for both cases. The DAL cases are either three or
five times coarser than the DNS case in all three spatial directions. Similar
to the forced isotropic turbulence case, we utilize a down-sampling method
where the velocity field at a certain point of the LES is the same as the
(interpolated) DNS velocity field at the same point. The details of the
simulations are given in Table 1.

Case | Method | $Re_{\tau}$ | $\Delta_{x}^{+}$ | $\min(\Delta_{y}^{+})$ | $\max(\Delta_{y}^{+})$ | $\Delta_{z}^{+}$
---|---|---|---|---|---|---
CH180DNS | DNS | $186$ | 6.1 | 0.19 | 8.2 | 6.1
CH180DAL3 | DAL | $186$ | 18.3 | 0.61 | 24.3 | 18.3
CH550DNS | DNS | $547$ | 6.7 | 0.21 | 9.6 | 3.4
CH550DAL3 | DAL | $547$ | 20.2 | 0.66 | 28.7 | 10.1
CH550DAL5 | DAL | $547$ | 33.6 | 0.97 | 49.5 | 16.8
CH550DSM3 | DSM | $547$ | 20.2 | 0.66 | 28.7 | 10.1
CH550AMD3 | AMD | $547$ | 20.2 | 0.66 | 28.7 | 10.1
CH550NM3 | NM | $547$ | 20.2 | 0.66 | 28.7 | 10.1
Table 1: Tabulated list of cases for the channel flow case. Cases are named as
CH[$Re_{\tau}$][Method][Coarsening Factor], where method is either DNS or DAL
and coarsening factor is the ratio of grid points from DNS to DAL in one
spatial direction.

(a)  (b)

(c)  (d)

Figure 3: Time evolution of the (a) integrated numerical error
$\mathcal{E}^{\text{num}}$ and wall-normal distance dependent numerical error
$\mathcal{E}^{\text{num}}_{y}$ at (b) $y^{+}\approx 15$ and (c) $y/h\approx 1$
for CH180DAL3 ( ), CH550DAL3 ( ), and CH550DAL5 ( ). Vertical lines ( ) divide
the initial response region, growth region, and saturation region. (d)
$\mathcal{E}^{\text{num}}_{y}$ normalized with the mean velocity profile as a
function of wall-normal height at $tu_{\tau}^{2}/\nu\approx 1$ for CH550DAL3.

In order to evaluate the numerical error, we use the same approach as in
section 3.1 and supply the SGS stress term using DAL for the three different
cases. Figure 3(a) shows the evolution of the integrated numerical error as a
function of time. When the error is normalized with
$\Delta_{\text{DNS}}/\Delta_{\text{DAL}}$, we see that the error evolution of
the three cases collapses. The linear dependence of the error on grid
resolution and its independence on Reynolds number is similar to the LES error
scalings found in Lozano-Durán and Bae [21]. This study did not distinguish
between modeling and numerical errors and considered errors in resulting
statistics. However, this shows that either the numerical errors are dominant
in the total LES error or the modeling error is proportional to the numerical
error, demonstrating that numerical errors are a significant source of error
in LES.

For all three cases, after the initial linear growth, the errors assume a
polynomial growth with a growth rate of 0.8, until saturation. This growth
rate is larger than that observed in the forced isotropic turbulence case. The
faster growth in error can be attributed to the non-uniform mesh in the wall-
normal direction and, incidentally, the modifications in the Poisson solver,
which contribute towards additional sources of numerical error.

Figures 3(b) and (c) show $\mathcal{E}^{\text{num}}_{y}$ for a given time at
different wall-normal heights, $y^{+}\approx 15$ and $y/\delta=1$. The linear
dependence on the factor $\Delta_{\text{DNS}}/\Delta_{\text{DAL}}$ is true for
all wall-normal heights, similar to the integrated error. The growth rates are
similar for the two wall-normal heights, with the error accumulating slightly
slower in the outer region, compared to the error at $y^{+}\approx 15$. This
shows that the numerical error will accumulate faster closer to the wall.
Additionally, the errors are weakly dependent on Reynolds number close to the
wall, but become independent of Reynolds number away from the wall. Thus, the
relative error close to the wall will increase with Reynolds number but stay
constant with Reynolds number in the outer region. Moreover, the relative
$y$-dependent numerical error at each wall-normal height is given in figure
3(d). While only one time step is shown, all temporal locations show similar
trends as a function of wall-normal distance. This shows that the relative
numerical error is the largest close to the wall and decreases with distance
from the wall. This large and fast accumulation of errors in the near-wall
region is known to be problematic for SGS models and in particular for wall
models for LES.

Finally, we plot the modeling and total errors for the DSM and AMD model as
well as the NM case in figure 4. The errors show that the AMD model has the
largest modeling error $\mathcal{E}_{mod}$, followed by the DSM and the NM
case. However, the total error $\mathcal{E}_{tot}$ shows the smallest error
for the AMD case, followed by the DSM and NM. This demonstrates that the
modeling error and the numerical error cancel each other, as widely accepted
in the literature.

(a)  (b)

Figure 4: (a) $\mathcal{E}^{\text{mod}}_{y}$ and (b)
$\mathcal{E}^{\text{tot}}_{y}$ normalized with the mean velocity profile as a
function of wall-normal height at $tu_{\tau}^{2}/\nu\approx 1$ for CH550DSM3 (
), CH550AMD3 ( ), and CH550NM3 ( ).

##  4 Conclusions

Large-eddy simulation has become a viable tool in performing simulations of
turbulent flow for industrial and engineering applications. In order for LES
to be regarded as a useful tool, the errors associated with it must be
carefully studied. In particular, the errors associated with the numerical
method and the choice of SGS model must be clearly distinguished and evaluated
at grid resolutions relevant to practical LES applications. Moreover, for
wall-bounded flows, it is important to evaluate the error as a function of
wall-normal distance to account for near-wall effects.

In this paper, we introduce a formulation that allows the accurate assessment
of numerical errors at grid resolutions typical of LES. DNS-aided LES is used
to compute exact SGS stress terms in the absence of numerical error; that is,
by running a DNS and LES with the equivalent initial condition, the DNS flow
field at each time step is used to generate the exact SGS stress term for the
corresponding LES flow field. This requires an explicit filter operator which
is often lacking in traditional LES. DAL allows the evaluation of numerical
errors as a function of time. Moreover, it allows the distinction of numerical
and modeling errors for LES, which is often not straightforward.

Using this method, we compute the numerical errors of two cases: forced
isotropic turbulence and turbulent channel flow. The evaluation of numerical
errors of the Fourier discretization shows that the errors are machine
precision zero, as expected. This also verifies that DAL indeed produces zero
modeling error. In the case of numerical error for the finite difference case,
the error follows a polynomial growth until saturation. The growth rate is
indicative of how the LES solution is correlated to the corresponding DNS
solution. The results verify that the numerical errors are largest closer to
the wall. Furthermore, it also shows that numerical errors and modeling errors
partially counteract each other. The findings imply that SGS models are not
expected to perform well in all numerical methods and will perform better when
coupled with numerical methods where the numerical errors are counterbalanced
with the modeling errors. In particular, supervised learning SGS models
trained on filtered DNS data do not account for numerical errors, and thus do
not perform as expected in _a posteriori_ testing. The latter observation also
has implications for wall modeling in LES, in which wall models become
sensitive to the flow details close to the wall where the numerical errors are
the largest.

The method introduced here can be helpful in two main directions. Firstly, the
method can inform the development of new SGS models in conjunction with a
target numerical method such that the effect of counterbalancing the numerical
and modeling errors is maximized. Since the method here provides the numerical
errors as a function of space and time, this information can be used to
develop SGS models specifically designed to cancel numerical errors. Secondly,
the method can similarly aid in the development and assessment of wall models,
which remains a pacing item to achieve practical LES.

## Acknowledgments

H.J.B acknowledges support from the National Science Foundation under grant
No.2152705.

##  Appendix A Commutation of filter operator for finite difference

For the filter operator to commute with the differentiation operator
numerically, we need, in matrix form,

$D^{\text{LES}}F\vec{u}=FD^{\text{DNS}}\vec{u},$ (10)

where a nonzero matrix $F\in\mathbb{R}^{N\times mN}$ is the filter matrix,
$D^{\text{DNS}}\in\mathbb{R}^{mN\times mN}$ and
$D^{\text{LES}}\in\mathbb{R}^{N\times N}$ are the one-dimensional
differentiation matrix for DNS and LES respectively. The factor $m$ denotes
the ratio between the DNS and LES grid sizes, and for practical cases, $m>1$.
For simplicity, we consider the central finite difference in a uniform grid,
which gives

$D^{\text{DNS}}_{i,j}=\begin{cases}1/\Delta,&i=j+1\\\ -1/\Delta,&i=j-1\\\
0,&\text{otherwise}\end{cases},\quad
D^{\text{LES}}_{i,j}=\begin{cases}1/m\Delta,&i=j+1\\\ -1/m\Delta,&i=j-1\\\
0,&\text{otherwise}\end{cases},$ (11)

where $\Delta$ is the grid size for DNS.

Solving Eq. (10) as $F=D^{\text{LES}-1}FD^{\text{DNS}}$, we get

$F_{i,j}=\begin{cases}\displaystyle{m\sum_{k=(i+1)/2}^{N/2}\left(-F_{2k,j-1}+F_{2k,j+1}\right)},&i\text{
odd}\\\
\displaystyle{m\sum_{k=1/2}^{(i-1)/2}\left(F_{2k,j-1}-F_{2k,j+1}\right)},&i\text{
even}\end{cases}$ (12)

We can rewrite this equation as a linear equation $M\vec{f}=\vec{f}$ by
transforming matrix $F$ as a vector $\vec{f}\in\mathbb{R}^{mN^{2}}$, such that
$\vec{f}_{mN(i-1)+j}=F_{i,j}$. Thus, there exists a filter that commutes with
the finite difference if and only if $M$ has an eigenvalue of unity. While
there may be a specific choice of $m$ and $N$ such that this is possible (e.g.
$m=1$), for an arbitrary grid, this is not always true.

## References

  * Smagorinsky [1963] J. Smagorinsky,  General circulation experiments with the primitive equations,  Mon. Weather Rev. 91 (1963) 99–164.
  * Schumann [1975] U. Schumann,  Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli,  J. Comput. Phys. 18 (1975) 376–404.
  * Kraichnan [1976] R. H. Kraichnan,  Eddy viscosity in two and three dimensions,  J. Atmos. Sci. 33 (1976) 1521–1536.
  * Métais and Lesieur [1992] O. Métais, M. Lesieur,  Spectral large-eddy simulation of isotropic and stably stratified turbulence,  J. Fluid Mech. 239 (1992) 157–194.
  * Rozema et al. [2015] W. Rozema, H. J. Bae, P. Moin, R. Verstappen,  Minimum-dissipation models for large-eddy simulation,  Phys. Fluids 27 (2015) 085107.
  * Yoshizawa [1982] A. Yoshizawa,  A statistically-derived subgrid model for the large-eddy simulation of turbulence,  Phys. Fluids 25 (1982) 1532–1538.
  * Leith [1990] C. E. Leith,  Stochastic backscatter in a subgrid-scale model: Plane shear mixing layer,  Phys. Fluids A 2 (1990) 297–299.
  * Clark et al. [1979] R. A. Clark, J. H. Ferziger, W. C. Reynolds,  Evaluation of subgrid-scale models using an accurately simulated turbulent flow,  J. Fluid Mech. 91 (1979) 1–16.
  * Leonard [1997] A. Leonard,  Large-eddy simulation of chaotic convection and beyond,  in: 35th Aerospace Sciences Meeting and Exhibit, 1997, p. 204.
  * Geurts [1997] B. J. Geurts,  Inverse modeling for large-eddy simulation,  Phys. Fluids 9 (1997) 3585–3587.
  * Zang et al. [1993] Y. Zang, R. L. Street, J. R. Koseff,  A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows,  Phys. Fluids A 5 (1993) 3186–3196.
  * Vreman et al. [1994] B. Vreman, B. Geurts, H. Kuerten,  On the formulation of the dynamic mixed subgrid-scale model,  Phys. Fluids 6 (1994) 4057–4059.
  * Germano et al. [1991] M. Germano, U. Piomelli, P. Moin, W. H. Cabot,  A dynamic subgrid-scale eddy viscosity model,  Phys. Fluids A 3 (1991) 1760.
  * Germano [1992] M. Germano,  Turbulence: the filtering approach,  J. Fluid Mech. 238 (1992) 325–336.
  * Lilly [1992] D. K. Lilly,  A proposed modification of the Germano subgrid-scale closure method,  Phys. Fluids A 4 (1992) 633–635.
  * Kobayashi and Shimomura [2001] H. Kobayashi, Y. Shimomura,  The performance of dynamic subgrid-scale models in the large eddy simulation of rotating homogeneous turbulence,  Phys. Fluids 13 (2001) 2350–2360.
  * Vreman et al. [1996] B. Vreman, B. Geurts, H. Kuerten,  Large-eddy simulation of the temporal mixing layer using the Clark model,  Theor. Comp. Fluid Dyn. 8 (1996) 309–324.
  * Vreman et al. [1997] B. Vreman, B. Geurts, H. Kuerten,  Large-eddy simulation of the turbulent mixing layer,  J. Fluid Mech. 339 (1997) 357–390.
  * Piomelli et al. [1988] U. Piomelli, P. Moin, J. H. Ferziger,  Model consistency in large eddy simulation of turbulent channel flows,  Phys. Fluids 31 (1988) 1884–1891.
  * Chung and McKeon [2010] D. Chung, B. J. McKeon,  Large-eddy simulation of large-scale structures in long channel flow,  J. Fluid Mech. 661 (2010) 341–364.
  * Lozano-Durán and Bae [2019] A. Lozano-Durán, H. J. Bae,  Error scaling of large-eddy simulation in the outer region of wall-bounded turbulence,  J. Comput. Phys. 392 (2019) 532–555.
  * Bonnet et al. [1998] J. Bonnet, R. Moser, W. Rodi,  A selection of test cases for the validation of large eddy simulations of turbulent flows,  AGARD advisory report 345 (1998) 1–35.
  * Vreman et al. [1996] B. Vreman, B. Geurts, H. Kuerten,  Comparision of numerical schemes in large-eddy simulation of the temporal mixing layer,  Int. J. Numer. Methods Fluids 22 (1996) 297–311.
  * Geurts and Fröhlich [2002] B. J. Geurts, J. Fröhlich,  A framework for predicting accuracy limitations in large-eddy simulation,  Phys. Fluids 14 (2002) L41–L44.
  * Meyers et al. [2003] J. Meyers, B. J. Geurts, M. Baelmans,  Database analysis of errors in large-eddy simulation,  Phys. Fluids 15 (2003) 2740–2755.
  * Meyers and Sagaut [2007] J. Meyers, P. Sagaut,  Is plane-channel flow a friendly case for the testing of large-eddy simulation subgrid-scale models?,  Phys. Fluids 19 (2007) 048105.
  * Meyers [2011] J. Meyers,  Error-landscape assessment of large-eddy simulations: A review of the methodology,  J. Sci. Comput. 49 (2011) 65–77.
  * Trias et al. [2017] F. X. Trias, A. Gorobets, M. H. Silvis, R. W. C. P. Verstappen, A. Oliva,  A new subgrid characteristic length for turbulence simulations on anisotropic grids,  Phys. Fluids 29 (2017) 115109.
  * Winckelmans et al. [2001] G. S. Winckelmans, A. A. Wray, O. V. Vasilyev, H. Jeanmart,  Explicit-filtering large-eddy simulation using the tensor-diffusivity model supplemented by a dynamic smagorinsky term,  Phys. Fluids 13 (2001) 1385–1403.
  * Stolz et al. [2001] S. Stolz, N. A. Adams, L. Kleiser,  An approximate deconvolution model for large-eddy simulation with application to incompressible wall-bounded flows,  Phys. Fluids 13 (2001) 997–1015.
  * Vasilyev et al. [1998] O. V. Vasilyev, T. S. Lund, P. Moin,  A general class of commutative filters for LES in complex geometries,  J. Comput. Phys. 146 (1998) 82–104.
  * Lund [2003] T. S. Lund,  The use of explicit filters in large eddy simulation,  Comput. Math. App. 46 (2003) 603–616.
  * Gullbrand and Chow [2003] J. Gullbrand, F. K. Chow,  The effect of numerical errors and turbulence models in large-eddy simulations of channel flow, with and without explicit filtering,  J. Fluid Mech. 495 (2003) 323–341.
  * Bose et al. [2010] S. T. Bose, P. Moin, D. You,  Grid-independent large-eddy simulation using explicit filtering,  Phys. Fluids 22 (2010) 105103.
  * Bae and Lozano-Durán [2017] H. J. Bae, A. Lozano-Durán,  Towards exact subgrid-scale models for explicitly filtered large-eddy simulation of wall-bounded flows,  Center for Turbulence Research - Annual Research Briefs (2017) 207–214.
  * Bae and Lozano-Durán [2018] H. J. Bae, A. Lozano-Durán,  DNS-aided explicitly filtered LES of channel flow,  Center for Turbulence Research - Annual Research Briefs (2018) 197–207.
  * Kutz [2017] J. N. Kutz,  Deep learning in fluid dynamics,  J. Fluid Mech. 814 (2017) 1–4.
  * Brenner et al. [2019] M. P. Brenner, J. D. Eldredge, J. B. Freund,  Perspective on machine learning for advancing fluid mechanics,  Phys. Rev. Fluids 4 (2019) 100501.
  * Duraisamy et al. [2019] K. Duraisamy, G. Iaccarino, H. Xiao,  Turbulence modeling in the age of data,  Annu. Rev. Fluid Mech. 51 (2019) 357–377.
  * Brunton et al. [2020] S. L. Brunton, B. R. Noack, P. Koumoutsakos,  Machine learning for fluid mechanics,  Annu. Rev. Fluid Mech. 52 (2020) 477–508.
  * Maulik et al. [2019] R. Maulik, O. San, A. Rasheed, P. Vedula,  Subgrid modelling for two-dimensional turbulence using neural networks,  J. Fluid Mech. 858 (2019) 122–144.
  * Guan et al. [2022] Y. Guan, A. Chattopadhyay, A. Subel, P. Hassanzadeh,  Stable a posteriori LES of 2D turbulence using convolutional neural networks: Backscattering analysis and generalization to higher $Re$ via transfer learning,  J. Comput. Phys. 458 (2022) 111090.
  * Zhou et al. [2019] Z. Zhou, G. He, S. Wang, G. Jin,  Subgrid-scale model for large-eddy simulation of isotropic turbulent flows using an artificial neural network,  Comput. Fluids 195 (2019) 104319.
  * Xie et al. [2020] C. Xie, J. Wang, E. Weinan,  Modeling subgrid-scale forces by spatial artificial neural networks in large eddy simulation of turbulence,  Phys. Rev. Fluids 5 (2020) 054606.
  * Frezat et al. [2021] H. Frezat, G. Balarac, J. Le Sommer, R. Fablet, R. Lguensat,  Physical invariance in neural networks for subgrid-scale scalar flux modeling,  Phys. Rev. Fluids 6 (2021) 024607.
  * Novati et al. [2021] G. Novati, H. L. de Laroussilhe, P. Koumoutsakos,  Automating turbulence modelling by multi-agent reinforcement learning,  Nat. Mach. Intell. 3 (2021) 87–96.
  * Wang et al. [2018] Z. Wang, K. Luo, D. Li, J. Tan, J. Fan,  Investigations of data-driven closure for subgrid-scale stress in large-eddy simulation,  Phys. Fluids 30 (2018) 125101.
  * Beck et al. [2019] A. Beck, D. Flad, C.-D. Munz,  Deep neural networks for data-driven les closure models,  J. Comput. Phys. 398 (2019) 108910.
  * Gamahara and Hattori [2017] M. Gamahara, Y. Hattori,  Searching for turbulence models by artificial neural network,  Phys. Rev. Fluids 2 (2017) 054604.
  * Park and Choi [2021] J. Park, H. Choi,  Toward neural-network-based large eddy simulation: Application to turbulent channel flow,  J. Fluid Mech. 914 (2021).
  * Bae and Koumoutsakos [2022] H. J. Bae, P. Koumoutsakos,  Scientific multi-agent reinforcement learning for wall-models of turbulent flows,  Nat. Commun. 13 (2022) 1–9.
  * Lozano-Durán and Bae [2020] A. Lozano-Durán, H. J. Bae,  Self-critical machine-learning wall-modeled LES for external aerodynamics,  arXiv preprint arXiv:2012.10005 (2020).
  * Nadiga and Livescu [2007] B. T. Nadiga, D. Livescu,  Instability of the perfect subgrid model in implicit-filtering large eddy simulation of geostrophic turbulence,  Phys. Rev. E 75 (2007) 046303.
  * Marsden et al. [2002] A. L. Marsden, O. V. Vasilyev, P. Moin,  Construction of commutative filters for les on unstructured meshes,  J Comput. Phys. 175 (2002) 584–603.
  * Lundgren [2003] T. S. Lundgren,  Linearly forced isotropic turbulence,  Center for Turbulence Research - Annual Research Briefs (2003) 461–473.
  * Nastac et al. [2017] G. Nastac, J. W. Labahn, L. Magri, M. Ihme,  Lyapunov exponent as a metric for assessing the dynamic content and predictability of large-eddy simulations,  Phys. Rev. Fluids 2 (2017) 094606.
  * Germano [1986] M. Germano,  Differential filters for the large eddy numerical simulation of turbulent flows,  Phys. Fluids 29 (1986) 1755–1757.
  * Orlandi [2000] P. Orlandi, Fluid Flow Phenomena: A Numerical Toolkit, 1, Springer, 2000\.
  * Kim and Moin [1985] J. Kim, P. Moin,  Application of a fractional-step method to incompressible Navier-Stokes equations,  J. Comput. Phys. 59 (1985) 308–323.
  * Wray [1990] A. A. Wray, Minimal-storage time advancement schemes for spectral methods, Technical Report, NASA Ames Research Center, 1990.
  * Lozano-Durán and Bae [2016] A. Lozano-Durán, H. J. Bae,  Turbulent channel with slip boundaries as a benchmark for subgrid-scale models in LES,  Center for Turbulence Research - Annual Research Briefs (2016) 97–103.
  * Bae et al. [2018] H. J. Bae, A. Lozano-Duran, S. Bose, P. Moin,  Turbulence intensities in large-eddy simulation of wall-bounded flows,  Phys. Rev. Fluids 3 (2018) 014610.
  * Bae et al. [2019] H. J. Bae, A. Lozano-Durán, S. T. Bose, P. Moin,  Dynamic slip wall model for large-eddy simulation,  J. Fluid Mech. 859 (2019) 400–432.
  * Lozano-Durán and Jiménez [2014] A. Lozano-Durán, J. Jiménez,  Effect of the computational domain on direct simulations of turbulent channels up to $Re_{\tau}$ = 4200,  Phys. Fluids 26 (2014) 011702.