Truncation Error-Based Anisotropic p-Adaptation for Unsteady Flows for High-Order Discontinuous Galerkin Methods

2025-05-06 0 0 2.26MB 23 页 10玖币
侵权投诉
Truncation Error-Based Anisotropic p-Adaptation for Unsteady Flows for
High-Order Discontinuous Galerkin Methods
Andr´es M. Rueda-Ram´ıreza,, Gerasimos Ntoukasb, Gonzalo Rubiob,c,, Eusebio Valerob,c, Esteban
Ferrerb,c
aDepartment of Mathematics and Computer Science, University of Cologne, 50931 Cologne, Germany
bETSIAE-UPM (School of Aeronautics), Universidad Polit´ecnica de Madrid, Plaza de Cardenal Cisneros 3, 28040 Madrid,
Spain
cCenter for Computational Simulation, Universidad Polit´ecnica de Madrid, Campus de Montegancedo, Boadilla del Monte,
28660 Madrid, Spain
Abstract
In this work, we extend the τ-estimation method to unsteady problems and use it to adapt the polynomial
degree for high-order discontinuous Galerkin simulations of unsteady flows. The adaptation is local and
anisotropic and allows capturing relevant unsteady flow features while enhancing the accuracy of time evolv-
ing functionals (e.g., lift, drag). To achieve an efficient and unsteady truncation error-based p-adaptation
scheme, we first revisit the definition of the truncation error, studying the effect of the treatment of the
mass matrix arising from the temporal term. Secondly, we extend the τ-estimation strategy to unsteady
problems. Finally, we present and compare two adaptation strategies for unsteady problems: the dynamic
and static p-adaptation methods. In the first one (dynamic) the error is measured periodically during a
simulation and the polynomial degree is adapted immediately after every estimation procedure. In the
second one (static) the error is also measured periodically, but only one p-adaptation process is performed
after several estimation stages, using a combination of the periodic error measures. The static p-adaptation
strategy is suitable for time-periodic flows, while the dynamic one can be generalized to any flow evolution.
We consider two test cases to evaluate the efficiency of the proposed p-adaptation strategies. The first
one considers the compressible Euler equations to simulate the advection of a density pulse. The second
one solves the compressible Navier-Stokes equations to simulate the flow around a cylinder at Re=100. The
local and anisotropic adaptation enables significant reductions in the number of degrees of freedom with
respect to uniform refinement, leading to speed-ups of up to ×4.5 for the Euler test case and ×2.2 for the
Navier-Stokes test case.
Keywords: High-order discontinuous Galerkin, Anisotropic p-adaptation, Unsteady p-adaptation,
Compressible flows.
1. Introduction
High-order DG methods are expected to be the engine of the next generation of CFD codes [1]. Com-
monly, high-order DG schemes are formulated as multi-domain spectral methods. As a result, besides the
increased accuracy of spectral (high-order) methods, they also provide a compact stencil and therefore a
local character, a feature that makes them highly parallelizable and flexible for complex 3D geometries [2,3].
Moreover, DG methods can handle non-conforming meshes with hanging nodes and/or different polynomial
degrees efficiently [46], which makes them well suited for mesh adaptation strategies.
Among the different DG formulations, the discontinuous Galerkin spectral element method (DGSEM)
[7,8] is a nodal (collocation) version of the DG method that uses tensor-product Lagrange basis functions
Corresponding authors:
Email addresses: aruedara@uni-koeln.de (Andr´es M. Rueda-Ram´ırez ), g.rubio@upm.es (Gonzalo Rubio )
Preprint submitted to Journal of Applied Mathematics and Computation October 10, 2022
arXiv:2210.03523v1 [math.NA] 7 Oct 2022
(traditionally in hexahedra) and stores the data at the Gauss or Gauss-Lobatto nodes of the quadrature
rule. The use of a quadrature rule with the same number of nodes as the approximate solution equips the
DGSEM with a diagonal mass matrix and very cheap-to-compute operators. In fact, the computational
cost of the DGSEM has been estimated to be a factor of four smaller than other conventional DG methods
[9]. In addition, since the DGSEM uses tensor-product bases, it can handle p-anisotropic discretizations
efficiently [10,11], i.e. discretizations that have different polynomial degrees in each coordinate direction.
For all those properties, the DGSEM has been used in a wide range of applications, including the simulation
of incompressible Navier-Stokes, compressible Navier-Stokes, Cahn Hilliard equation, or multi-phase flows
(see [12] and references therein).
In their famous review paper, Wang et al. [3] point out that one of the challenges the high-order commu-
nity must address to impact the design process and replace traditional low-order codes is the development
of efficient mesh adaptation strategies. The idea behind these strategies is to reduce the number of degrees
of freedom (DOFs) while maintaining high accuracy, which translates into shorter computational times and
reduced storage requirements. Local adaptation can be performed by subdividing or merging elements (h-
adaptation), by enriching or reducing the polynomial degree in certain elements (p-adaptation), by relocating
the position of the nodes in a mesh (r-adaptation). For all these strategies it is of paramount importance
to identify the flow regions that require refinement or coarsening with a local error estimation.
Adaptation strategies have been classified according to the type of error measure that is employed as
feature-based adaptation, adjoint-based adaptation, and local error-based adaptation. A comparison of these
three approaches was performed by Fraysse et al. [13] for finite volume approximations and by Kompenhans
et al. [14] and Naddei et al. [15] for high-order DG methods. The feature-based adaptation is the classical
approach and uses easy-to-compute error measures that depend on the flow features. They rely on the
assumption that high errors are expected where the flow is more difficult to resolve. Hence, refinement is
predicted where high velocity, density or pressure gradients are identified [16,17]. For DG discretizations,
an easy-to-compute feature-based adaptation criterion is the assessment of jumps across element interfaces
[1820]. The main disadvantage of these methods is that there is no direct relation between the adaptation
criterion and the numerical errors and thus the accuracy is not easily predictable. Additionally, the only
way to solve steady-state problems is to adapt iteratively.
A second and more sophisticated approach is known as adjoint-based adaptation. In this approach, a
functional target is defined (e.g. drag or lift in external flow aerodynamics) and the adjoint problem is solved
to obtain a spatial distribution of the functional error, which is then used to adapt the mesh. This technique
was originally developed for structural analysis using FEM by Babuˇska and Miller [21,22], and has been
used recently for adaptation strategies in DG methods [2325]. The main drawback of this approach is the
high computational cost to solve the adjoint problem and the storage requirements needed to save the error
estimators, especially in unsteady flows. Moreover, only the error of the functional analyzed is guaranteed
to be reduced, whereas the error of other functionals may deteriorate.
A computationally more efficient alternative is the local error-based adaptation, which is based on the
assessment of any measurable (not feature-based) local error in all the cells of the domain [24]. The local
error-based adaptation methods are interesting since, in contrast to feature-based methods, they provide a
way to predict and control the overall accuracy, and are computationally cheaper than adjoint-based schemes
[10,14]. For those reasons, local error-based adaptation strategies are retained in this work. A large amount
of effort has been invested in the development of reliable local error-based adaptation methods. Estimations
of the local discretization error have been used by Mavriplis [26,27] to develop hp-adaptation techniques for
the spectral element method. Residual-based p-adaptation is also a local error-based adaptation method,
which uses the residual to measure how accurate is the local approximation. This method was originally
developed for Finite Elements (FE) and has been successfully used with DG methods [15,23]. In the case
of modal (hierarchical) DG methods, a possibility is to employ low cost error estimates that take advantage
of the modal approximation to drive p-adaptation procedures, such as the Variational Multiscale (VMS)
indicator by Kuru and De la Llave Plata [28], or the spectral decay indicator by Persson and Peraire [17].
In this work, we favor truncation error estimators, another local error-based alternative to drive a mesh
adaptation method. The truncation error is related to the discretization error through the Discretization
Error Transport Equation [29], where it acts as a local source term. This relation makes it useful as an
2
indicator for mesh adaptation methods [30,31], since refining the mesh where the truncation error is high
reduces the discretization error in all the mesh [32], with an additional advantage: truncation error estimation
requires less computational effort than adjoint methods. Finally, it has been shown that controlling the
truncation error targets the numerical accuracy of all functionals at once [10,33], ensuring that adapting a
mesh using the truncation error leads necessarily to an error decrease in any other functional (e.g. lift or
drag).
The τ-estimation method proposed by Brandt [34], which estimates the local truncation error by injecting
a fine grid solution into coarser meshes, has been used to perform local error-based mesh adaptation in low-
order schemes [13,31,3538]. Rubio et al. [39] extended the τ-estimation approach to high-order methods
using a continuous Chebyshev collocation method. Later, Rubio et al. [32] applied it to DGSEM discretiza-
tions. Kompenhans et al. [10] applied the τ-estimation approach to perform steady-state p-adaptation using
the Euler and Navier-Stokes equations, and showed that a reduction of the truncation error increases the
numerical accuracy of all functionals at once. Furthermore, Kompenhans et al. [14] also showed that trun-
cation error-based adaptation can exhibit better performance than feature-based adaptation. In contrast
to most local error-based adaptation methods, where multiple error estimation and adaptation stages are
needed in a steady-state solution, the τ-estimation method generates a unique prediction of what polynomial
degree is needed for a desired truncation error threshold. Therefore, in steady-state the adaptation strategy
is to converge a high-order approximation (reference mesh) to a specified global residual and then to perform
a single error estimation followed by a corresponding p-adaptation process. Besides, the truncation error is
known to decay exponentially in smooth solutions [10,32]. Therefore, if the estimation is good, it is possible
to extrapolate the behavior and predict the polynomial degree needed for a desired error threshold [10,40].
Being a relatively recent technique, p-adaptation methods that use τ-estimators have only been applied
to steady-state solutions with high-order methods [10,11]. In this work, we propose a methodology to
extend this technique to unsteady problems. The rest of this work is organized as follows. First, in Section
2.2, we show that the truncation error can be formulated in several forms, depending on the choice of
the continuous and discrete partial differential operators. Herein, a thorough analysis of the formulation
that is traditionally used in the DG community [10,14,32,40,41] is presented, a new formulation for the
truncation error is proposed, and a comparative analysis of both techniques is detailed. Second, in Section 2.3
we provide two strategies to estimate the truncation error in unsteady problems, both of which are derived
from the τ-estimation method. The first strategy is directly derived from the variational DG formulation
while the second strategy uses a dual time-stepping pseudo-time discretization. Third, in Section 2.4, two p-
adaptation algorithms for unsteady problems are proposed, which use an unsteady τ-estimation method: (i)
a dynamic p-adaptation strategy, which performs several stages of estimation and p-adaptation throughout
a simulation and second, and (ii) a static p-adaptation strategy, which performs several truncation error
estimation stages, but only one p-adaptation stage. Finally, the methods are applied to unsteady problems
modeled by the compressible Euler and Navier-Stokes equations in Section 3, and a detailed analysis of their
performance is presented. The most important findings of this paper are summarized in Section 4.
2. Numerical Methods
2.1. The Discontinuous Galerkin Spectral Element Method
We consider the approximation of systems of conservation laws,
tq+F(q)=0,in Ω,(1)
subject to appropriate boundary conditions, where qis the state vector of conserved variables, and F(q)=
fis the continuous partial differential operator, where
fis a flux block vector, which depends on q.
In an advection-diffusion conservation law, such as the Navier-Stokes equations, the flux vector can be
written as
f=
fa(q)
fν(q,
q),(2)
3
where
fais the advective flux and
fνis the diffusive flux. Because of the dependency of the diffusive flux
on
q, (1) is a second order PDE. Following Arnold et al. [42], (1) can be rewritten as a first-order system,
tq+
fa(q)
fν(q,
g)=0,in Ω,
q=
g,in Ω.
(3a)
(3b)
To obtain the DGSEM-version of (3), the computational domain is subdivided into non-overlapping
hexahedral elements, all variables are approximated by piece-wise Lagrange interpolating polynomials of
degree Nthat are continuous in each element, but allowed to be discontinuous across element interfaces:
qqN,
f
fNand
g
gN. Furthermore, (3a) and (3b) are multiplied by an arbitrary polynomial
(test function) of degree N, the derivative terms are integrated by parts, and all integrals are evaluated
numerically with a quadrature rule of N+1 points, to obtain
JjwjtqN
jN
e
fN
φjdΩe+N
e
ˆ
fφjdSe=0,
N
eqN
φjdΩe+N
eφjˆ
q
ndSe=Jjwj
gN
j
(4a)
(4b)
for each degree of freedom of each element. In (4), ˆ
fand ˆ
qare the numerical traces of the flux and the solu-
tion, respectively, the functions φjare the so-called basis functions, which are tensor product expansions of
the Lagrange interpolating polynomials, the Jjare the Jacobians of the geometry transformation with which
the mesh is created, and the wjare the weights of the quadrature rule. The derivation of (4) is given in [8,43].
The discretization of the system can be compactly written as
MNdQN
dt+F
F
FN(QN)=0,(5)
where MNis the mass matrix of the system, QNa vector with all the unknowns, and F
F
FNthe discrete partial
differential operator. The mass matrix of the DGSEM is diagonal so, (5) is often rewritten as
dQN
dt+(MN)1F
F
FN(QN)=0.(6)
2.2. Formulation of the Truncation Error
The truncation error is defined as the difference between the discrete partial differential operator and
the continuous partial differential operator, both applied to the exact solution of the problem. This is often
known as Generalized Truncation Error Expresion (GTEE) [44,45]. For the problem at hand, we take
the difference between the discrete equation (6) applied to the sampled continuous solution, INq, and the
sampled continuous equation (1) to obtain:
dINq
dtINtq+(MN)1F
F
FN(INq)INF(q)=˜
τ
τ
τN,(7)
where ˜
τ
τ
τNis the truncation error and INis a restriction/prolongation operator used to project the solution
from one space to another. Here we are simply sampling the continuous solution into our discrete space.
Assuming that the restriction operator commutes with the time derivative (INtq=dINq
dt), the two first
terms in (7) cancel out,
˜
τ
τ
τN=(MN)1F
F
FN(INq)INF(q).(8)
The discretization error is defined as the difference between the approximate solution and the exact
solution to the problem, i.e.,
N=QNINq. Taking the difference between (6) and (1) sampled into the
discrete space, we get the following.
dQN
dtINtq+(MN)1F
F
FN(QN)INF(q)=0.(9)
4
By the linearity of the time derivative and using the definition of the discretization error,
d
N
dt=(MN)1F
F
FN(
N+INq)+INF(q).(10)
Now, for simplicity, we consider that F
F
FNis a linear operator (it can be linearized otherwise to achieve a
similar result; see, e.g., [46]),
d
N
dt+(MN)1F
F
FN(
N)=(MN)1F
F
FN(INq)+INF(q)=˜
τ
τ
τN.(11)
Therefore, for linear operators, the discretization error is governed by the same equation as the numerical
solution with the addition of the truncation error as a source term. This equation is known as Discrete
Error Transport Equation (DETE). As can be seen, the truncation error is interesting for mesh adaptation,
as it acts as a source for the generation of discretization error.
The truncation error definition, (8), can be re-scaled with the mass matrix:
τ
τ
τN=MN˜
τ
τ
τN=F
F
FN(INq)MNINF(q),(12)
which is equivalent to defining the truncation error as the projection of the difference between discrete and
continuous operators on the individual basis functions of the finite element subspace. With this definition,
the DETE reads:
MNd
N
dt+F
F
FN(
N)=F
F
FN(INq)+MNINF(q)=τ
τ
τN.(13)
The definition (12) has previously been used in [10,14,32,40,41], and will be called traditional formu-
lation in this work. The definition (8) has been recently used in [33] and will be called, in this work, new
formulation. The main advantage of the new formulation is that it is directly related to functional errors,
as shown in [33].
Note that previous works in the context of the DGSEM [10,14,32,33,40,41] focused primarily on
steady-state problems; therefore, the second term in the RHS of (8) and (12) was zero.
We have derived two formulations of the truncation error and showed that each leads to a different
version of the DETE: (11) and (13). Although both formulations appear to be similar, some remarks can
be made about their properties.
1. The traditional version of the truncation error (12) acts as a source term for the discretization error,
after projecting it point-wise on the basis functions φj, which build the finite element subspace, as
shown in (13). On the other hand, the new approximation of the truncation error (8) acts directly as
a source term of the pointwise values of the discretization error; see (11).
2. Since DGSEM is a collocation method and the mass matrix is a diagonal matrix containing the
mapping Jacobian and quadrature weights, see, for example, [12], each form of the truncation error
can be obtained from the other by scaling it point-wise with Jjwj, see (12). In other words, the main
difference between both formulations is the weight that they give to the element size.
3. Due to the strong similarities between the two truncation error approximations, the anisotropic prop-
erties and the possibility to estimate the error in a multigrid cycle (see [41]) hold for both error
measures.
Due to the similarities of both truncation error formulations, the tilde notation will be dropped in next
sections, and the expressions will hold for both, unless the contrary is explicitly stated.
On a final note, in [32] the concept of isolated truncation error was introduced. While the standard
non-isolated truncation error uses all terms appearing in the discrete discontinuous Galerkin variational
formulation, (4), the isolated truncation error replaces the surface numerical flux functions by simple evalu-
ations of the flux with the inner solution of each element. The isolated truncation error has some advantages
for mesh adaptation, as shown in [14].
5
摘要:

TruncationError-BasedAnisotropicp-AdaptationforUnsteadyFlowsforHigh-OrderDiscontinuousGalerkinMethodsAndresM.Rueda-Ramreza,‡,GerasimosNtoukasb,GonzaloRubiob,c,‡,EusebioValerob,c,EstebanFerrerb,caDepartmentofMathematicsandComputerScience,UniversityofCologne,50931Cologne,GermanybETSIAE-UPM(Schoolof...

展开>> 收起<<
Truncation Error-Based Anisotropic p-Adaptation for Unsteady Flows for High-Order Discontinuous Galerkin Methods.pdf

共23页,预览5页

还剩页未读, 继续阅读

声明:本站为文档C2C交易模式,即用户上传的文档直接被用户下载,本站只是中间服务平台,本站所有文档下载所得的收益归上传人(含作者)所有。玖贝云文库仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对上载内容本身不做任何修改或编辑。若文档所含内容侵犯了您的版权或隐私,请立即通知玖贝云文库,我们立即给予删除!
分类:图书资源 价格:10玖币 属性:23 页 大小:2.26MB 格式:PDF 时间:2025-05-06

开通VIP享超值会员特权

  • 多端同步记录
  • 高速下载文档
  • 免费文档工具
  • 分享文档赚钱
  • 每日登录抽奖
  • 优质衍生服务
/ 23
客服
关注